1 | // Copyright (C) 2002, International Business Machines |
---|
2 | // Corporation and others. All Rights Reserved. |
---|
3 | |
---|
4 | #ifdef NDEBUG |
---|
5 | #undef NDEBUG |
---|
6 | #endif |
---|
7 | |
---|
8 | #include "ClpConfig.h" |
---|
9 | #include "CoinPragma.hpp" |
---|
10 | #include <cassert> |
---|
11 | #include <cstdio> |
---|
12 | #include <cstdlib> |
---|
13 | #include <cmath> |
---|
14 | #include <cfloat> |
---|
15 | #include <string> |
---|
16 | #include <iostream> |
---|
17 | |
---|
18 | #include "CoinMpsIO.hpp" |
---|
19 | #include "CoinPackedMatrix.hpp" |
---|
20 | #include "CoinPackedVector.hpp" |
---|
21 | #include "CoinStructuredModel.hpp" |
---|
22 | #include "CoinHelperFunctions.hpp" |
---|
23 | #include "CoinTime.hpp" |
---|
24 | |
---|
25 | #include "ClpFactorization.hpp" |
---|
26 | #include "ClpSimplex.hpp" |
---|
27 | #include "ClpSimplexOther.hpp" |
---|
28 | #include "ClpSimplexNonlinear.hpp" |
---|
29 | #include "ClpInterior.hpp" |
---|
30 | #include "ClpLinearObjective.hpp" |
---|
31 | #include "ClpDualRowSteepest.hpp" |
---|
32 | #include "ClpDualRowDantzig.hpp" |
---|
33 | #include "ClpPrimalColumnSteepest.hpp" |
---|
34 | #include "ClpPrimalColumnDantzig.hpp" |
---|
35 | #include "ClpParameters.hpp" |
---|
36 | #include "ClpNetworkMatrix.hpp" |
---|
37 | #include "ClpPlusMinusOneMatrix.hpp" |
---|
38 | #include "MyMessageHandler.hpp" |
---|
39 | #include "MyEventHandler.hpp" |
---|
40 | |
---|
41 | #include "ClpPresolve.hpp" |
---|
42 | #include "Idiot.hpp" |
---|
43 | |
---|
44 | |
---|
45 | //############################################################################# |
---|
46 | |
---|
47 | |
---|
48 | // Function Prototypes. Function definitions is in this file. |
---|
49 | void testingMessage( const char * const msg ); |
---|
50 | #if UFL_BARRIER |
---|
51 | static int barrierAvailable=1; |
---|
52 | static std::string nameBarrier="barrier-UFL"; |
---|
53 | #elif WSSMP_BARRIER |
---|
54 | static int barrierAvailable=2; |
---|
55 | static std::string nameBarrier="barrier-WSSMP"; |
---|
56 | #elif MUMPS_BARRIER |
---|
57 | static int barrierAvailable=3; |
---|
58 | static std::string nameBarrier="barrier-MUMPS"; |
---|
59 | #else |
---|
60 | static int barrierAvailable=0; |
---|
61 | static std::string nameBarrier="barrier-slow"; |
---|
62 | #endif |
---|
63 | #define NUMBER_ALGORITHMS 12 |
---|
64 | // If you just want a subset then set some to 1 |
---|
65 | static int switchOff[NUMBER_ALGORITHMS]={0,0,0,0,0,0,0,0,0,0,0,0}; |
---|
66 | // shortName - 0 no , 1 yes |
---|
67 | ClpSolve setupForSolve(int algorithm, std::string & nameAlgorithm, |
---|
68 | int shortName) |
---|
69 | { |
---|
70 | ClpSolve solveOptions; |
---|
71 | /* algorithms are |
---|
72 | 0 barrier |
---|
73 | 1 dual with volumne crash |
---|
74 | 2,3 dual with and without crash |
---|
75 | 4,5 primal with and without |
---|
76 | 6,7 automatic with and without |
---|
77 | 8,9 primal with idiot 1 and 5 |
---|
78 | 10,11 primal with 70, dual with volume |
---|
79 | */ |
---|
80 | switch (algorithm) { |
---|
81 | case 0: |
---|
82 | if (shortName) |
---|
83 | nameAlgorithm="ba"; |
---|
84 | else |
---|
85 | nameAlgorithm="nameBarrier"; |
---|
86 | solveOptions.setSolveType(ClpSolve::useBarrier); |
---|
87 | if (barrierAvailable==1) |
---|
88 | solveOptions.setSpecialOption(4,4); |
---|
89 | else if (barrierAvailable==2) |
---|
90 | solveOptions.setSpecialOption(4,2); |
---|
91 | break; |
---|
92 | case 1: |
---|
93 | #ifdef COIN_HAS_VOL |
---|
94 | if (shortName) |
---|
95 | nameAlgorithm="du-vol-50"; |
---|
96 | else |
---|
97 | nameAlgorithm="dual-volume-50"; |
---|
98 | solveOptions.setSolveType(ClpSolve::useDual); |
---|
99 | solveOptions.setSpecialOption(0,2,50); // volume |
---|
100 | #else |
---|
101 | solveOptions.setSolveType(ClpSolve::notImplemented); |
---|
102 | #endif |
---|
103 | break; |
---|
104 | case 2: |
---|
105 | if (shortName) |
---|
106 | nameAlgorithm="du-cr"; |
---|
107 | else |
---|
108 | nameAlgorithm="dual-crash"; |
---|
109 | solveOptions.setSolveType(ClpSolve::useDual); |
---|
110 | solveOptions.setSpecialOption(0,1); |
---|
111 | break; |
---|
112 | case 3: |
---|
113 | if (shortName) |
---|
114 | nameAlgorithm="du"; |
---|
115 | else |
---|
116 | nameAlgorithm="dual"; |
---|
117 | solveOptions.setSolveType(ClpSolve::useDual); |
---|
118 | break; |
---|
119 | case 4: |
---|
120 | if (shortName) |
---|
121 | nameAlgorithm="pr-cr"; |
---|
122 | else |
---|
123 | nameAlgorithm="primal-crash"; |
---|
124 | solveOptions.setSolveType(ClpSolve::usePrimal); |
---|
125 | solveOptions.setSpecialOption(1,1); |
---|
126 | break; |
---|
127 | case 5: |
---|
128 | if (shortName) |
---|
129 | nameAlgorithm="pr"; |
---|
130 | else |
---|
131 | nameAlgorithm="primal"; |
---|
132 | solveOptions.setSolveType(ClpSolve::usePrimal); |
---|
133 | break; |
---|
134 | case 6: |
---|
135 | if (shortName) |
---|
136 | nameAlgorithm="au-cr"; |
---|
137 | else |
---|
138 | nameAlgorithm="either-crash"; |
---|
139 | solveOptions.setSolveType(ClpSolve::automatic); |
---|
140 | solveOptions.setSpecialOption(1,1); |
---|
141 | break; |
---|
142 | case 7: |
---|
143 | if (shortName) |
---|
144 | nameAlgorithm="au"; |
---|
145 | else |
---|
146 | nameAlgorithm="either"; |
---|
147 | solveOptions.setSolveType(ClpSolve::automatic); |
---|
148 | break; |
---|
149 | case 8: |
---|
150 | if (shortName) |
---|
151 | nameAlgorithm="pr-id-1"; |
---|
152 | else |
---|
153 | nameAlgorithm="primal-idiot-1"; |
---|
154 | solveOptions.setSolveType(ClpSolve::usePrimalorSprint); |
---|
155 | solveOptions.setSpecialOption(1,2,1); // idiot |
---|
156 | break; |
---|
157 | case 9: |
---|
158 | if (shortName) |
---|
159 | nameAlgorithm="pr-id-5"; |
---|
160 | else |
---|
161 | nameAlgorithm="primal-idiot-5"; |
---|
162 | solveOptions.setSolveType(ClpSolve::usePrimalorSprint); |
---|
163 | solveOptions.setSpecialOption(1,2,5); // idiot |
---|
164 | break; |
---|
165 | case 10: |
---|
166 | if (shortName) |
---|
167 | nameAlgorithm="pr-id-70"; |
---|
168 | else |
---|
169 | nameAlgorithm="primal-idiot-70"; |
---|
170 | solveOptions.setSolveType(ClpSolve::usePrimalorSprint); |
---|
171 | solveOptions.setSpecialOption(1,2,70); // idiot |
---|
172 | break; |
---|
173 | case 11: |
---|
174 | #ifdef COIN_HAS_VOL |
---|
175 | if (shortName) |
---|
176 | nameAlgorithm="du-vol"; |
---|
177 | else |
---|
178 | nameAlgorithm="dual-volume"; |
---|
179 | solveOptions.setSolveType(ClpSolve::useDual); |
---|
180 | solveOptions.setSpecialOption(0,2,3000); // volume |
---|
181 | #else |
---|
182 | solveOptions.setSolveType(ClpSolve::notImplemented); |
---|
183 | #endif |
---|
184 | break; |
---|
185 | default: |
---|
186 | abort(); |
---|
187 | } |
---|
188 | if (shortName) { |
---|
189 | // can switch off |
---|
190 | if (switchOff[algorithm]) |
---|
191 | solveOptions.setSolveType(ClpSolve::notImplemented); |
---|
192 | } |
---|
193 | return solveOptions; |
---|
194 | } |
---|
195 | static void printSol(ClpSimplex & model) |
---|
196 | { |
---|
197 | int numberRows = model.numberRows(); |
---|
198 | int numberColumns = model.numberColumns(); |
---|
199 | |
---|
200 | double * rowPrimal = model.primalRowSolution(); |
---|
201 | double * rowDual = model.dualRowSolution(); |
---|
202 | double * rowLower = model.rowLower(); |
---|
203 | double * rowUpper = model.rowUpper(); |
---|
204 | |
---|
205 | int iRow; |
---|
206 | double objValue = model.getObjValue(); |
---|
207 | printf("Objvalue %g Rows (%d)\n",objValue,numberRows); |
---|
208 | for (iRow=0;iRow<numberRows;iRow++) { |
---|
209 | printf("%d primal %g dual %g low %g up %g\n", |
---|
210 | iRow,rowPrimal[iRow],rowDual[iRow], |
---|
211 | rowLower[iRow],rowUpper[iRow]); |
---|
212 | } |
---|
213 | double * columnPrimal = model.primalColumnSolution(); |
---|
214 | double * columnDual = model.dualColumnSolution(); |
---|
215 | double * columnLower = model.columnLower(); |
---|
216 | double * columnUpper = model.columnUpper(); |
---|
217 | double offset; |
---|
218 | //const double * gradient = model.objectiveAsObject()->gradient(&model, |
---|
219 | // columnPrimal,offset,true,1); |
---|
220 | const double * gradient = model.objective(columnPrimal,offset); |
---|
221 | int iColumn; |
---|
222 | objValue = -offset-model.objectiveOffset(); |
---|
223 | printf("offset %g (%g)\n",offset,model.objectiveOffset()); |
---|
224 | printf("Columns (%d)\n",numberColumns); |
---|
225 | for (iColumn=0;iColumn<numberColumns;iColumn++) { |
---|
226 | printf("%d primal %g dual %g low %g up %g\n", |
---|
227 | iColumn,columnPrimal[iColumn],columnDual[iColumn], |
---|
228 | columnLower[iColumn],columnUpper[iColumn]); |
---|
229 | objValue += columnPrimal[iColumn]*gradient[iColumn]; |
---|
230 | if (fabs(columnPrimal[iColumn]*gradient[iColumn])>1.0e-8) |
---|
231 | printf("obj -> %g gradient %g\n",objValue,gradient[iColumn]); |
---|
232 | } |
---|
233 | printf("Computed objective %g\n",objValue); |
---|
234 | } |
---|
235 | |
---|
236 | void usage(const std::string& key) |
---|
237 | { |
---|
238 | std::cerr |
---|
239 | <<"Undefined parameter \"" <<key <<"\".\n" |
---|
240 | <<"Correct usage: \n" |
---|
241 | <<" clp [-dirSample=V1] [-dirNetlib=V2] [-netlib]\n" |
---|
242 | <<" where:\n" |
---|
243 | <<" -dirSample: directory containing mps test files\n" |
---|
244 | <<" Default value V1=\"../../Data/Sample\"\n" |
---|
245 | <<" -dirNetlib: directory containing netlib files\"\n" |
---|
246 | <<" Default value V2=\"../../Data/Netlib\"\n" |
---|
247 | <<" -netlib\n" |
---|
248 | <<" If specified, then netlib testset run as well as the nitTest.\n"; |
---|
249 | } |
---|
250 | |
---|
251 | //---------------------------------------------------------------- |
---|
252 | int mainTest (int argc, const char *argv[],int algorithm, |
---|
253 | ClpSimplex empty, bool doPresolve, int switchOffValue,bool doVector) |
---|
254 | { |
---|
255 | int i; |
---|
256 | |
---|
257 | if (switchOffValue>0) { |
---|
258 | // switch off some |
---|
259 | int iTest; |
---|
260 | for (iTest=0;iTest<NUMBER_ALGORITHMS;iTest++) { |
---|
261 | #ifndef NDEBUG |
---|
262 | int bottom = switchOffValue%10; |
---|
263 | assert (bottom==0||bottom==1); |
---|
264 | #endif |
---|
265 | switchOffValue/= 10; |
---|
266 | switchOff[iTest]=0; |
---|
267 | } |
---|
268 | } |
---|
269 | |
---|
270 | // define valid parameter keywords |
---|
271 | std::set<std::string> definedKeyWords; |
---|
272 | definedKeyWords.insert("-dirSample"); |
---|
273 | definedKeyWords.insert("-dirNetlib"); |
---|
274 | definedKeyWords.insert("-netlib"); |
---|
275 | |
---|
276 | // Create a map of parameter keys and associated data |
---|
277 | std::map<std::string,std::string> parms; |
---|
278 | for ( i=1; i<argc; i++ ) { |
---|
279 | std::string parm(argv[i]); |
---|
280 | std::string key,value; |
---|
281 | std::string::size_type eqPos = parm.find('='); |
---|
282 | |
---|
283 | // Does parm contain and '=' |
---|
284 | if ( eqPos==std::string::npos ) { |
---|
285 | //Parm does not contain '=' |
---|
286 | key = parm; |
---|
287 | } |
---|
288 | else { |
---|
289 | key=parm.substr(0,eqPos); |
---|
290 | value=parm.substr(eqPos+1); |
---|
291 | } |
---|
292 | |
---|
293 | // Is specifed key valid? |
---|
294 | if ( definedKeyWords.find(key) == definedKeyWords.end() ) { |
---|
295 | // invalid key word. |
---|
296 | // Write help text |
---|
297 | usage(key); |
---|
298 | return 1; |
---|
299 | } |
---|
300 | parms[key]=value; |
---|
301 | } |
---|
302 | |
---|
303 | const char dirsep = CoinFindDirSeparator(); |
---|
304 | // Set directory containing mps data files. |
---|
305 | std::string dirSample; |
---|
306 | if (parms.find("-dirSample") != parms.end()) |
---|
307 | dirSample=parms["-dirSample"]; |
---|
308 | else |
---|
309 | dirSample = dirsep == '/' ? "../../Data/Sample/" : "..\\..\\Data\\Sample\\"; |
---|
310 | |
---|
311 | // Set directory containing netlib data files. |
---|
312 | std::string dirNetlib; |
---|
313 | if (parms.find("-dirNetlib") != parms.end()) |
---|
314 | dirNetlib=parms["-dirNetlib"]; |
---|
315 | else |
---|
316 | dirNetlib = dirsep == '/' ? "../../Data/Netlib/" : "..\\..\\Data\\Netlib\\"; |
---|
317 | if (!empty.numberRows()) { |
---|
318 | testingMessage( "Testing ClpSimplex\n" ); |
---|
319 | ClpSimplexUnitTest(dirSample); |
---|
320 | } |
---|
321 | if (parms.find("-netlib") != parms.end()||empty.numberRows()) |
---|
322 | { |
---|
323 | unsigned int m; |
---|
324 | |
---|
325 | // Define test problems: |
---|
326 | // mps names, |
---|
327 | // maximization or minimization, |
---|
328 | // Number of rows and columns in problem, and |
---|
329 | // objective function value |
---|
330 | std::vector<std::string> mpsName; |
---|
331 | std::vector<bool> min; |
---|
332 | std::vector<int> nRows; |
---|
333 | std::vector<int> nCols; |
---|
334 | std::vector<double> objValue; |
---|
335 | std::vector<double> objValueTol; |
---|
336 | // 100 added means no presolve |
---|
337 | std::vector<int> bestStrategy; |
---|
338 | if(empty.numberRows()) { |
---|
339 | std::string alg; |
---|
340 | for (int iTest=0;iTest<NUMBER_ALGORITHMS;iTest++) { |
---|
341 | ClpSolve solveOptions=setupForSolve(iTest,alg,0); |
---|
342 | printf("%d %s ",iTest,alg.c_str()); |
---|
343 | if (switchOff[iTest]) |
---|
344 | printf("skipped by user\n"); |
---|
345 | else if(solveOptions.getSolveType()==ClpSolve::notImplemented) |
---|
346 | printf("skipped as not available\n"); |
---|
347 | else |
---|
348 | printf("will be tested\n"); |
---|
349 | } |
---|
350 | } |
---|
351 | if (!empty.numberRows()) { |
---|
352 | mpsName.push_back("25fv47"); |
---|
353 | min.push_back(true); |
---|
354 | nRows.push_back(822); |
---|
355 | nCols.push_back(1571); |
---|
356 | objValueTol.push_back(1.E-10); |
---|
357 | objValue.push_back(5.5018458883E+03); |
---|
358 | bestStrategy.push_back(0); |
---|
359 | mpsName.push_back("80bau3b");min.push_back(true);nRows.push_back(2263);nCols.push_back(9799);objValueTol.push_back(1.e-10);objValue.push_back(9.8722419241E+05);bestStrategy.push_back(3); |
---|
360 | mpsName.push_back("blend");min.push_back(true);nRows.push_back(75);nCols.push_back(83);objValueTol.push_back(1.e-10);objValue.push_back(-3.0812149846e+01);bestStrategy.push_back(3); |
---|
361 | mpsName.push_back("pilotnov");min.push_back(true);nRows.push_back(976);nCols.push_back(2172);objValueTol.push_back(1.e-10);objValue.push_back(-4.4972761882e+03);bestStrategy.push_back(3); |
---|
362 | mpsName.push_back("maros-r7");min.push_back(true);nRows.push_back(3137);nCols.push_back(9408);objValueTol.push_back(1.e-10);objValue.push_back(1.4971851665e+06);bestStrategy.push_back(2); |
---|
363 | mpsName.push_back("pilot");min.push_back(true);nRows.push_back(1442);nCols.push_back(3652);objValueTol.push_back(1.e-5);objValue.push_back(/*-5.5740430007e+02*/-557.48972927292);bestStrategy.push_back(3); |
---|
364 | mpsName.push_back("pilot4");min.push_back(true);nRows.push_back(411);nCols.push_back(1000);objValueTol.push_back(5.e-5);objValue.push_back(-2.5811392641e+03);bestStrategy.push_back(3); |
---|
365 | mpsName.push_back("pilot87");min.push_back(true);nRows.push_back(2031);nCols.push_back(4883);objValueTol.push_back(1.e-4);objValue.push_back(3.0171072827e+02);bestStrategy.push_back(0); |
---|
366 | mpsName.push_back("adlittle");min.push_back(true);nRows.push_back(57);nCols.push_back(97);objValueTol.push_back(1.e-10);objValue.push_back(2.2549496316e+05);bestStrategy.push_back(3); |
---|
367 | mpsName.push_back("afiro");min.push_back(true);nRows.push_back(28);nCols.push_back(32);objValueTol.push_back(1.e-10);objValue.push_back(-4.6475314286e+02);bestStrategy.push_back(3); |
---|
368 | mpsName.push_back("agg");min.push_back(true);nRows.push_back(489);nCols.push_back(163);objValueTol.push_back(1.e-10);objValue.push_back(-3.5991767287e+07);bestStrategy.push_back(3); |
---|
369 | mpsName.push_back("agg2");min.push_back(true);nRows.push_back(517);nCols.push_back(302);objValueTol.push_back(1.e-10);objValue.push_back(-2.0239252356e+07);bestStrategy.push_back(3); |
---|
370 | mpsName.push_back("agg3");min.push_back(true);nRows.push_back(517);nCols.push_back(302);objValueTol.push_back(1.e-10);objValue.push_back(1.0312115935e+07);bestStrategy.push_back(4); |
---|
371 | mpsName.push_back("bandm");min.push_back(true);nRows.push_back(306);nCols.push_back(472);objValueTol.push_back(1.e-10);objValue.push_back(-1.5862801845e+02);bestStrategy.push_back(2); |
---|
372 | mpsName.push_back("beaconfd");min.push_back(true);nRows.push_back(174);nCols.push_back(262);objValueTol.push_back(1.e-10);objValue.push_back(3.3592485807e+04);bestStrategy.push_back(0); |
---|
373 | mpsName.push_back("bnl1");min.push_back(true);nRows.push_back(644);nCols.push_back(1175);objValueTol.push_back(1.e-10);objValue.push_back(1.9776295615E+03);bestStrategy.push_back(3); |
---|
374 | mpsName.push_back("bnl2");min.push_back(true);nRows.push_back(2325);nCols.push_back(3489);objValueTol.push_back(1.e-10);objValue.push_back(1.8112365404e+03);bestStrategy.push_back(3); |
---|
375 | mpsName.push_back("boeing1");min.push_back(true);nRows.push_back(/*351*/352);nCols.push_back(384);objValueTol.push_back(1.e-10);objValue.push_back(-3.3521356751e+02);bestStrategy.push_back(3); |
---|
376 | mpsName.push_back("boeing2");min.push_back(true);nRows.push_back(167);nCols.push_back(143);objValueTol.push_back(1.e-10);objValue.push_back(-3.1501872802e+02);bestStrategy.push_back(3); |
---|
377 | mpsName.push_back("bore3d");min.push_back(true);nRows.push_back(234);nCols.push_back(315);objValueTol.push_back(1.e-10);objValue.push_back(1.3730803942e+03);bestStrategy.push_back(3); |
---|
378 | mpsName.push_back("brandy");min.push_back(true);nRows.push_back(221);nCols.push_back(249);objValueTol.push_back(1.e-10);objValue.push_back(1.5185098965e+03);bestStrategy.push_back(3); |
---|
379 | mpsName.push_back("capri");min.push_back(true);nRows.push_back(272);nCols.push_back(353);objValueTol.push_back(1.e-10);objValue.push_back(2.6900129138e+03);bestStrategy.push_back(3); |
---|
380 | mpsName.push_back("cycle");min.push_back(true);nRows.push_back(1904);nCols.push_back(2857);objValueTol.push_back(1.e-9);objValue.push_back(-5.2263930249e+00);bestStrategy.push_back(3); |
---|
381 | mpsName.push_back("czprob");min.push_back(true);nRows.push_back(930);nCols.push_back(3523);objValueTol.push_back(1.e-10);objValue.push_back(2.1851966989e+06);bestStrategy.push_back(3); |
---|
382 | mpsName.push_back("d2q06c");min.push_back(true);nRows.push_back(2172);nCols.push_back(5167);objValueTol.push_back(1.e-7);objValue.push_back(122784.21557456);bestStrategy.push_back(0); |
---|
383 | mpsName.push_back("d6cube");min.push_back(true);nRows.push_back(416);nCols.push_back(6184);objValueTol.push_back(1.e-7);objValue.push_back(3.1549166667e+02);bestStrategy.push_back(3); |
---|
384 | mpsName.push_back("degen2");min.push_back(true);nRows.push_back(445);nCols.push_back(534);objValueTol.push_back(1.e-10);objValue.push_back(-1.4351780000e+03);bestStrategy.push_back(3); |
---|
385 | mpsName.push_back("degen3");min.push_back(true);nRows.push_back(1504);nCols.push_back(1818);objValueTol.push_back(1.e-10);objValue.push_back(-9.8729400000e+02);bestStrategy.push_back(2); |
---|
386 | mpsName.push_back("dfl001");min.push_back(true);nRows.push_back(6072);nCols.push_back(12230);objValueTol.push_back(1.e-5);objValue.push_back(1.1266396047E+07);bestStrategy.push_back(5); |
---|
387 | mpsName.push_back("e226");min.push_back(true);nRows.push_back(224);nCols.push_back(282);objValueTol.push_back(1.e-10);objValue.push_back(-1.8751929066e+01+7.113);bestStrategy.push_back(3); // The correct answer includes -7.113 term. This is a constant in the objective function. See line 1683 of the mps file. |
---|
388 | mpsName.push_back("etamacro");min.push_back(true);nRows.push_back(401);nCols.push_back(688);objValueTol.push_back(1.e-6);objValue.push_back(-7.5571521774e+02 );bestStrategy.push_back(3); |
---|
389 | mpsName.push_back("fffff800");min.push_back(true);nRows.push_back(525);nCols.push_back(854);objValueTol.push_back(1.e-6);objValue.push_back(5.5567961165e+05);bestStrategy.push_back(3); |
---|
390 | mpsName.push_back("finnis");min.push_back(true);nRows.push_back(498);nCols.push_back(614);objValueTol.push_back(1.e-6);objValue.push_back(1.7279096547e+05);bestStrategy.push_back(3); |
---|
391 | mpsName.push_back("fit1d");min.push_back(true);nRows.push_back(25);nCols.push_back(1026);objValueTol.push_back(1.e-10);objValue.push_back(-9.1463780924e+03);bestStrategy.push_back(3+100); |
---|
392 | mpsName.push_back("fit1p");min.push_back(true);nRows.push_back(628);nCols.push_back(1677);objValueTol.push_back(1.e-10);objValue.push_back(9.1463780924e+03);bestStrategy.push_back(5+100); |
---|
393 | mpsName.push_back("fit2d");min.push_back(true);nRows.push_back(26);nCols.push_back(10500);objValueTol.push_back(1.e-10);objValue.push_back(-6.8464293294e+04);bestStrategy.push_back(3+100); |
---|
394 | mpsName.push_back("fit2p");min.push_back(true);nRows.push_back(3001);nCols.push_back(13525);objValueTol.push_back(1.e-9);objValue.push_back(6.8464293232e+04);bestStrategy.push_back(5+100); |
---|
395 | mpsName.push_back("forplan");min.push_back(true);nRows.push_back(162);nCols.push_back(421);objValueTol.push_back(1.e-6);objValue.push_back(-6.6421873953e+02);bestStrategy.push_back(3); |
---|
396 | mpsName.push_back("ganges");min.push_back(true);nRows.push_back(1310);nCols.push_back(1681);objValueTol.push_back(1.e-5);objValue.push_back(-1.0958636356e+05);bestStrategy.push_back(3); |
---|
397 | mpsName.push_back("gfrd-pnc");min.push_back(true);nRows.push_back(617);nCols.push_back(1092);objValueTol.push_back(1.e-10);objValue.push_back(6.9022359995e+06);bestStrategy.push_back(3); |
---|
398 | mpsName.push_back("greenbea");min.push_back(true);nRows.push_back(2393);nCols.push_back(5405);objValueTol.push_back(1.e-10);objValue.push_back(/*-7.2462405908e+07*/-72555248.129846);bestStrategy.push_back(3); |
---|
399 | mpsName.push_back("greenbeb");min.push_back(true);nRows.push_back(2393);nCols.push_back(5405);objValueTol.push_back(1.e-10);objValue.push_back(/*-4.3021476065e+06*/-4302260.2612066);bestStrategy.push_back(3); |
---|
400 | mpsName.push_back("grow15");min.push_back(true);nRows.push_back(301);nCols.push_back(645);objValueTol.push_back(1.e-10);objValue.push_back(-1.0687094129e+08);bestStrategy.push_back(4+100); |
---|
401 | mpsName.push_back("grow22");min.push_back(true);nRows.push_back(441);nCols.push_back(946);objValueTol.push_back(1.e-10);objValue.push_back(-1.6083433648e+08);bestStrategy.push_back(4+100); |
---|
402 | mpsName.push_back("grow7");min.push_back(true);nRows.push_back(141);nCols.push_back(301);objValueTol.push_back(1.e-10);objValue.push_back(-4.7787811815e+07);bestStrategy.push_back(4+100); |
---|
403 | mpsName.push_back("israel");min.push_back(true);nRows.push_back(175);nCols.push_back(142);objValueTol.push_back(1.e-10);objValue.push_back(-8.9664482186e+05);bestStrategy.push_back(2); |
---|
404 | mpsName.push_back("kb2");min.push_back(true);nRows.push_back(44);nCols.push_back(41);objValueTol.push_back(1.e-10);objValue.push_back(-1.7499001299e+03);bestStrategy.push_back(3); |
---|
405 | mpsName.push_back("lotfi");min.push_back(true);nRows.push_back(154);nCols.push_back(308);objValueTol.push_back(1.e-10);objValue.push_back(-2.5264706062e+01);bestStrategy.push_back(3); |
---|
406 | mpsName.push_back("maros");min.push_back(true);nRows.push_back(847);nCols.push_back(1443);objValueTol.push_back(1.e-10);objValue.push_back(-5.8063743701e+04);bestStrategy.push_back(3); |
---|
407 | mpsName.push_back("modszk1");min.push_back(true);nRows.push_back(688);nCols.push_back(1620);objValueTol.push_back(1.e-10);objValue.push_back(3.2061972906e+02);bestStrategy.push_back(3); |
---|
408 | mpsName.push_back("nesm");min.push_back(true);nRows.push_back(663);nCols.push_back(2923);objValueTol.push_back(1.e-5);objValue.push_back(1.4076073035e+07);bestStrategy.push_back(2); |
---|
409 | mpsName.push_back("perold");min.push_back(true);nRows.push_back(626);nCols.push_back(1376);objValueTol.push_back(1.e-6);objValue.push_back(-9.3807580773e+03);bestStrategy.push_back(3); |
---|
410 | //mpsName.push_back("qap12");min.push_back(true);nRows.push_back(3193);nCols.push_back(8856);objValueTol.push_back(1.e-6);objValue.push_back(5.2289435056e+02);bestStrategy.push_back(3); |
---|
411 | //mpsName.push_back("qap15");min.push_back(true);nRows.push_back(6331);nCols.push_back(22275);objValueTol.push_back(1.e-10);objValue.push_back(1.0409940410e+03);bestStrategy.push_back(3); |
---|
412 | mpsName.push_back("recipe");min.push_back(true);nRows.push_back(92);nCols.push_back(180);objValueTol.push_back(1.e-10);objValue.push_back(-2.6661600000e+02);bestStrategy.push_back(3); |
---|
413 | mpsName.push_back("sc105");min.push_back(true);nRows.push_back(106);nCols.push_back(103);objValueTol.push_back(1.e-10);objValue.push_back(-5.2202061212e+01);bestStrategy.push_back(3); |
---|
414 | mpsName.push_back("sc205");min.push_back(true);nRows.push_back(206);nCols.push_back(203);objValueTol.push_back(1.e-10);objValue.push_back(-5.2202061212e+01);bestStrategy.push_back(3); |
---|
415 | mpsName.push_back("sc50a");min.push_back(true);nRows.push_back(51);nCols.push_back(48);objValueTol.push_back(1.e-10);objValue.push_back(-6.4575077059e+01);bestStrategy.push_back(3); |
---|
416 | mpsName.push_back("sc50b");min.push_back(true);nRows.push_back(51);nCols.push_back(48);objValueTol.push_back(1.e-10);objValue.push_back(-7.0000000000e+01);bestStrategy.push_back(3); |
---|
417 | mpsName.push_back("scagr25");min.push_back(true);nRows.push_back(472);nCols.push_back(500);objValueTol.push_back(1.e-10);objValue.push_back(-1.4753433061e+07);bestStrategy.push_back(3); |
---|
418 | mpsName.push_back("scagr7");min.push_back(true);nRows.push_back(130);nCols.push_back(140);objValueTol.push_back(1.e-6);objValue.push_back(-2.3313892548e+06);bestStrategy.push_back(3); |
---|
419 | mpsName.push_back("scfxm1");min.push_back(true);nRows.push_back(331);nCols.push_back(457);objValueTol.push_back(1.e-10);objValue.push_back(1.8416759028e+04);bestStrategy.push_back(3); |
---|
420 | mpsName.push_back("scfxm2");min.push_back(true);nRows.push_back(661);nCols.push_back(914);objValueTol.push_back(1.e-10);objValue.push_back(3.6660261565e+04);bestStrategy.push_back(3); |
---|
421 | mpsName.push_back("scfxm3");min.push_back(true);nRows.push_back(991);nCols.push_back(1371);objValueTol.push_back(1.e-10);objValue.push_back(5.4901254550e+04);bestStrategy.push_back(3); |
---|
422 | mpsName.push_back("scorpion");min.push_back(true);nRows.push_back(389);nCols.push_back(358);objValueTol.push_back(1.e-10);objValue.push_back(1.8781248227e+03);bestStrategy.push_back(3); |
---|
423 | mpsName.push_back("scrs8");min.push_back(true);nRows.push_back(491);nCols.push_back(1169);objValueTol.push_back(1.e-5);objValue.push_back(9.0429998619e+02);bestStrategy.push_back(2); |
---|
424 | mpsName.push_back("scsd1");min.push_back(true);nRows.push_back(78);nCols.push_back(760);objValueTol.push_back(1.e-10);objValue.push_back(8.6666666743e+00);bestStrategy.push_back(3+100); |
---|
425 | mpsName.push_back("scsd6");min.push_back(true);nRows.push_back(148);nCols.push_back(1350);objValueTol.push_back(1.e-10);objValue.push_back(5.0500000078e+01);bestStrategy.push_back(3+100); |
---|
426 | mpsName.push_back("scsd8");min.push_back(true);nRows.push_back(398);nCols.push_back(2750);objValueTol.push_back(1.e-10);objValue.push_back(9.0499999993e+02);bestStrategy.push_back(1+100); |
---|
427 | mpsName.push_back("sctap1");min.push_back(true);nRows.push_back(301);nCols.push_back(480);objValueTol.push_back(1.e-10);objValue.push_back(1.4122500000e+03);bestStrategy.push_back(3); |
---|
428 | mpsName.push_back("sctap2");min.push_back(true);nRows.push_back(1091);nCols.push_back(1880);objValueTol.push_back(1.e-10);objValue.push_back(1.7248071429e+03);bestStrategy.push_back(3); |
---|
429 | mpsName.push_back("sctap3");min.push_back(true);nRows.push_back(1481);nCols.push_back(2480);objValueTol.push_back(1.e-10);objValue.push_back(1.4240000000e+03);bestStrategy.push_back(3); |
---|
430 | mpsName.push_back("seba");min.push_back(true);nRows.push_back(516);nCols.push_back(1028);objValueTol.push_back(1.e-10);objValue.push_back(1.5711600000e+04);bestStrategy.push_back(3); |
---|
431 | mpsName.push_back("share1b");min.push_back(true);nRows.push_back(118);nCols.push_back(225);objValueTol.push_back(1.e-10);objValue.push_back(-7.6589318579e+04);bestStrategy.push_back(3); |
---|
432 | mpsName.push_back("share2b");min.push_back(true);nRows.push_back(97);nCols.push_back(79);objValueTol.push_back(1.e-10);objValue.push_back(-4.1573224074e+02);bestStrategy.push_back(3); |
---|
433 | mpsName.push_back("shell");min.push_back(true);nRows.push_back(537);nCols.push_back(1775);objValueTol.push_back(1.e-10);objValue.push_back(1.2088253460e+09);bestStrategy.push_back(3); |
---|
434 | mpsName.push_back("ship04l");min.push_back(true);nRows.push_back(403);nCols.push_back(2118);objValueTol.push_back(1.e-10);objValue.push_back(1.7933245380e+06);bestStrategy.push_back(3); |
---|
435 | mpsName.push_back("ship04s");min.push_back(true);nRows.push_back(403);nCols.push_back(1458);objValueTol.push_back(1.e-10);objValue.push_back(1.7987147004e+06);bestStrategy.push_back(3); |
---|
436 | mpsName.push_back("ship08l");min.push_back(true);nRows.push_back(779);nCols.push_back(4283);objValueTol.push_back(1.e-10);objValue.push_back(1.9090552114e+06);bestStrategy.push_back(3); |
---|
437 | mpsName.push_back("ship08s");min.push_back(true);nRows.push_back(779);nCols.push_back(2387);objValueTol.push_back(1.e-10);objValue.push_back(1.9200982105e+06);bestStrategy.push_back(2); |
---|
438 | mpsName.push_back("ship12l");min.push_back(true);nRows.push_back(1152);nCols.push_back(5427);objValueTol.push_back(1.e-10);objValue.push_back(1.4701879193e+06);bestStrategy.push_back(3); |
---|
439 | mpsName.push_back("ship12s");min.push_back(true);nRows.push_back(1152);nCols.push_back(2763);objValueTol.push_back(1.e-10);objValue.push_back(1.4892361344e+06);bestStrategy.push_back(2); |
---|
440 | mpsName.push_back("sierra");min.push_back(true);nRows.push_back(1228);nCols.push_back(2036);objValueTol.push_back(1.e-10);objValue.push_back(1.5394362184e+07);bestStrategy.push_back(3); |
---|
441 | mpsName.push_back("stair");min.push_back(true);nRows.push_back(357);nCols.push_back(467);objValueTol.push_back(1.e-10);objValue.push_back(-2.5126695119e+02);bestStrategy.push_back(3); |
---|
442 | mpsName.push_back("standata");min.push_back(true);nRows.push_back(360);nCols.push_back(1075);objValueTol.push_back(1.e-10);objValue.push_back(1.2576995000e+03);bestStrategy.push_back(3); |
---|
443 | //mpsName.push_back("standgub");min.push_back(true);nRows.push_back(362);nCols.push_back(1184);objValueTol.push_back(1.e-10);objValue.push_back(1257.6995); bestStrategy.push_back(3); |
---|
444 | mpsName.push_back("standmps");min.push_back(true);nRows.push_back(468);nCols.push_back(1075);objValueTol.push_back(1.e-10);objValue.push_back(1.4060175000E+03); bestStrategy.push_back(3); |
---|
445 | mpsName.push_back("stocfor1");min.push_back(true);nRows.push_back(118);nCols.push_back(111);objValueTol.push_back(1.e-10);objValue.push_back(-4.1131976219E+04);bestStrategy.push_back(3); |
---|
446 | mpsName.push_back("stocfor2");min.push_back(true);nRows.push_back(2158);nCols.push_back(2031);objValueTol.push_back(1.e-10);objValue.push_back(-3.9024408538e+04);bestStrategy.push_back(3); |
---|
447 | //mpsName.push_back("stocfor3");min.push_back(true);nRows.push_back(16676);nCols.push_back(15695);objValueTol.push_back(1.e-10);objValue.push_back(-3.9976661576e+04);bestStrategy.push_back(3); |
---|
448 | //mpsName.push_back("truss");min.push_back(true);nRows.push_back(1001);nCols.push_back(8806);objValueTol.push_back(1.e-10);objValue.push_back(4.5881584719e+05);bestStrategy.push_back(3); |
---|
449 | mpsName.push_back("tuff");min.push_back(true);nRows.push_back(334);nCols.push_back(587);objValueTol.push_back(1.e-10);objValue.push_back(2.9214776509e-01);bestStrategy.push_back(3); |
---|
450 | mpsName.push_back("vtpbase");min.push_back(true);nRows.push_back(199);nCols.push_back(203);objValueTol.push_back(1.e-10);objValue.push_back(1.2983146246e+05);bestStrategy.push_back(3); |
---|
451 | mpsName.push_back("wood1p");min.push_back(true);nRows.push_back(245);nCols.push_back(2594);objValueTol.push_back(5.e-5);objValue.push_back(1.4429024116e+00);bestStrategy.push_back(3); |
---|
452 | mpsName.push_back("woodw");min.push_back(true);nRows.push_back(1099);nCols.push_back(8405);objValueTol.push_back(1.e-10);objValue.push_back(1.3044763331E+00);bestStrategy.push_back(3); |
---|
453 | } else { |
---|
454 | // Just testing one |
---|
455 | mpsName.push_back(empty.problemName());min.push_back(true);nRows.push_back(-1); |
---|
456 | nCols.push_back(-1);objValueTol.push_back(1.e-10); |
---|
457 | objValue.push_back(0.0);bestStrategy.push_back(0); |
---|
458 | int iTest; |
---|
459 | std::string alg; |
---|
460 | for (iTest=0;iTest<NUMBER_ALGORITHMS;iTest++) { |
---|
461 | ClpSolve solveOptions=setupForSolve(iTest,alg,0); |
---|
462 | printf("%d %s ",iTest,alg.c_str()); |
---|
463 | if (switchOff[iTest]) |
---|
464 | printf("skipped by user\n"); |
---|
465 | else if(solveOptions.getSolveType()==ClpSolve::notImplemented) |
---|
466 | printf("skipped as not available\n"); |
---|
467 | else |
---|
468 | printf("will be tested\n"); |
---|
469 | } |
---|
470 | } |
---|
471 | |
---|
472 | double timeTaken =0.0; |
---|
473 | if( !barrierAvailable) |
---|
474 | switchOff[0]=1; |
---|
475 | // Loop once for each Mps File |
---|
476 | for (m=0; m<mpsName.size(); m++ ) { |
---|
477 | std::cerr <<" processing mps file: " <<mpsName[m] |
---|
478 | <<" (" <<m+1 <<" out of " <<mpsName.size() <<")" <<std::endl; |
---|
479 | |
---|
480 | ClpSimplex solutionBase=empty; |
---|
481 | std::string fn = dirNetlib+mpsName[m]; |
---|
482 | if (!empty.numberRows()||algorithm<6) { |
---|
483 | // Read data mps file, |
---|
484 | CoinMpsIO mps; |
---|
485 | int nerrors=mps.readMps(fn.c_str(),"mps"); |
---|
486 | if (nerrors) { |
---|
487 | std::cerr << "Error " << nerrors << " when reading model from " |
---|
488 | << fn.c_str() << "! Aborting tests.\n"; |
---|
489 | return 1; |
---|
490 | } |
---|
491 | solutionBase.loadProblem(*mps.getMatrixByCol(),mps.getColLower(), |
---|
492 | mps.getColUpper(), |
---|
493 | mps.getObjCoefficients(), |
---|
494 | mps.getRowLower(),mps.getRowUpper()); |
---|
495 | |
---|
496 | solutionBase.setDblParam(ClpObjOffset,mps.objectiveOffset()); |
---|
497 | } |
---|
498 | |
---|
499 | // Runs through strategies |
---|
500 | if (algorithm==6||algorithm==7) { |
---|
501 | // algorithms tested are at top of file |
---|
502 | double testTime[NUMBER_ALGORITHMS]; |
---|
503 | std::string alg[NUMBER_ALGORITHMS]; |
---|
504 | int iTest; |
---|
505 | for (iTest=0;iTest<NUMBER_ALGORITHMS;iTest++) { |
---|
506 | ClpSolve solveOptions=setupForSolve(iTest,alg[iTest],1); |
---|
507 | if (solveOptions.getSolveType()!=ClpSolve::notImplemented) { |
---|
508 | double time1 = CoinCpuTime(); |
---|
509 | ClpSimplex solution=solutionBase; |
---|
510 | if (solution.maximumSeconds()<0.0) |
---|
511 | solution.setMaximumSeconds(120.0); |
---|
512 | if (doVector) { |
---|
513 | ClpMatrixBase * matrix = solution.clpMatrix(); |
---|
514 | if (dynamic_cast< ClpPackedMatrix*>(matrix)) { |
---|
515 | ClpPackedMatrix * clpMatrix = dynamic_cast< ClpPackedMatrix*>(matrix); |
---|
516 | clpMatrix->makeSpecialColumnCopy(); |
---|
517 | } |
---|
518 | } |
---|
519 | solution.initialSolve(solveOptions); |
---|
520 | double time2 = CoinCpuTime()-time1; |
---|
521 | testTime[iTest]=time2; |
---|
522 | printf("Took %g seconds - status %d\n",time2,solution.problemStatus()); |
---|
523 | if (solution.problemStatus()) |
---|
524 | testTime[iTest]=1.0e20; |
---|
525 | } else { |
---|
526 | testTime[iTest]=1.0e30; |
---|
527 | } |
---|
528 | } |
---|
529 | int iBest=-1; |
---|
530 | double dBest=1.0e10; |
---|
531 | printf("%s",fn.c_str()); |
---|
532 | for (iTest=0;iTest<NUMBER_ALGORITHMS;iTest++) { |
---|
533 | if (testTime[iTest]<1.0e30) { |
---|
534 | printf(" %s %g", |
---|
535 | alg[iTest].c_str(),testTime[iTest]); |
---|
536 | if (testTime[iTest]<dBest) { |
---|
537 | dBest=testTime[iTest]; |
---|
538 | iBest=iTest; |
---|
539 | } |
---|
540 | } |
---|
541 | } |
---|
542 | printf("\n"); |
---|
543 | if (iBest>=0) |
---|
544 | printf("Best strategy for %s is %s (%d) which takes %g seconds\n", |
---|
545 | fn.c_str(),alg[iBest].c_str(),iBest,testTime[iBest]); |
---|
546 | else |
---|
547 | printf("No strategy finished in time limit\n"); |
---|
548 | continue; |
---|
549 | } |
---|
550 | double time1 = CoinCpuTime(); |
---|
551 | ClpSimplex solution=solutionBase; |
---|
552 | #if 0 |
---|
553 | solution.setOptimizationDirection(-1); |
---|
554 | { |
---|
555 | int j; |
---|
556 | double * obj = solution.objective(); |
---|
557 | int n=solution.numberColumns(); |
---|
558 | for (j=0;j<n;j++) |
---|
559 | obj[j] *= -1.0; |
---|
560 | } |
---|
561 | #endif |
---|
562 | ClpSolve::SolveType method; |
---|
563 | ClpSolve::PresolveType presolveType; |
---|
564 | ClpSolve solveOptions; |
---|
565 | if (doPresolve) |
---|
566 | presolveType=ClpSolve::presolveOn; |
---|
567 | else |
---|
568 | presolveType=ClpSolve::presolveOff; |
---|
569 | solveOptions.setPresolveType(presolveType,5); |
---|
570 | std::string nameAlgorithm; |
---|
571 | if (algorithm!=5) { |
---|
572 | if (algorithm==0) { |
---|
573 | method=ClpSolve::useDual; |
---|
574 | nameAlgorithm="dual"; |
---|
575 | } else if (algorithm==1) { |
---|
576 | method=ClpSolve::usePrimalorSprint; |
---|
577 | nameAlgorithm="primal"; |
---|
578 | } else if (algorithm==3) { |
---|
579 | method=ClpSolve::automatic; |
---|
580 | nameAlgorithm="either"; |
---|
581 | } else { |
---|
582 | nameAlgorithm="barrier-slow"; |
---|
583 | #ifdef UFL_BARRIER |
---|
584 | solveOptions.setSpecialOption(4,4); |
---|
585 | nameAlgorithm="barrier-UFL"; |
---|
586 | #endif |
---|
587 | #ifdef WSSMP_BARRIER |
---|
588 | solveOptions.setSpecialOption(4,2); |
---|
589 | nameAlgorithm="barrier-WSSMP"; |
---|
590 | #endif |
---|
591 | #ifdef MUMPS_BARRIER |
---|
592 | solveOptions.setSpecialOption(4,6); |
---|
593 | nameAlgorithm="barrier-MUMPS"; |
---|
594 | #endif |
---|
595 | method = ClpSolve::useBarrier; |
---|
596 | } |
---|
597 | solveOptions.setSolveType(method); |
---|
598 | } else { |
---|
599 | int iAlg = bestStrategy[m]; |
---|
600 | int presolveOff=iAlg/100; |
---|
601 | iAlg=iAlg % 100; |
---|
602 | if( !barrierAvailable&&iAlg==0) { |
---|
603 | if (nRows[m]!=2172) |
---|
604 | iAlg = 5; // try primal |
---|
605 | else |
---|
606 | iAlg=3; // d2q06c |
---|
607 | } |
---|
608 | solveOptions=setupForSolve(iAlg,nameAlgorithm,0); |
---|
609 | if (presolveOff) |
---|
610 | solveOptions.setPresolveType(ClpSolve::presolveOff); |
---|
611 | } |
---|
612 | if (doVector) { |
---|
613 | ClpMatrixBase * matrix = solution.clpMatrix(); |
---|
614 | if (dynamic_cast< ClpPackedMatrix*>(matrix)) { |
---|
615 | ClpPackedMatrix * clpMatrix = dynamic_cast< ClpPackedMatrix*>(matrix); |
---|
616 | clpMatrix->makeSpecialColumnCopy(); |
---|
617 | } |
---|
618 | } |
---|
619 | solution.initialSolve(solveOptions); |
---|
620 | double time2 = CoinCpuTime()-time1; |
---|
621 | timeTaken += time2; |
---|
622 | printf("%s took %g seconds using algorithm %s\n",fn.c_str(),time2,nameAlgorithm.c_str()); |
---|
623 | // Test objective solution value |
---|
624 | { |
---|
625 | double soln = solution.objectiveValue(); |
---|
626 | CoinRelFltEq eq(objValueTol[m]); |
---|
627 | std::cerr <<soln <<", " <<objValue[m] <<" diff "<< |
---|
628 | soln-objValue[m]<<std::endl; |
---|
629 | if(!eq(soln,objValue[m])) |
---|
630 | printf("** difference fails\n"); |
---|
631 | } |
---|
632 | } |
---|
633 | printf("Total time %g seconds\n",timeTaken); |
---|
634 | } |
---|
635 | else { |
---|
636 | testingMessage( "***Skipped Testing on netlib ***\n" ); |
---|
637 | testingMessage( "***use -netlib to test class***\n" ); |
---|
638 | } |
---|
639 | |
---|
640 | testingMessage( "All tests completed successfully\n" ); |
---|
641 | return 0; |
---|
642 | } |
---|
643 | |
---|
644 | |
---|
645 | // Display message on stdout and stderr |
---|
646 | void testingMessage( const char * const msg ) |
---|
647 | { |
---|
648 | std::cerr <<msg; |
---|
649 | //cout <<endl <<"*****************************************" |
---|
650 | // <<endl <<msg <<endl; |
---|
651 | } |
---|
652 | |
---|
653 | //-------------------------------------------------------------------------- |
---|
654 | // test factorization methods and simplex method and simple barrier |
---|
655 | void |
---|
656 | ClpSimplexUnitTest(const std::string & dirSample) |
---|
657 | { |
---|
658 | |
---|
659 | CoinRelFltEq eq(0.000001); |
---|
660 | |
---|
661 | { |
---|
662 | ClpSimplex solution; |
---|
663 | |
---|
664 | // matrix data |
---|
665 | //deliberate hiccup of 2 between 0 and 1 |
---|
666 | CoinBigIndex start[5]={0,4,7,8,9}; |
---|
667 | int length[5]={2,3,1,1,1}; |
---|
668 | int rows[11]={0,2,-1,-1,0,1,2,0,1,2}; |
---|
669 | double elements[11]={7.0,2.0,1.0e10,1.0e10,-2.0,1.0,-2.0,1,1,1}; |
---|
670 | CoinPackedMatrix matrix(true,3,5,8,elements,rows,start,length); |
---|
671 | |
---|
672 | // rim data |
---|
673 | double objective[7]={-4.0,1.0,0.0,0.0,0.0,0.0,0.0}; |
---|
674 | double rowLower[5]={14.0,3.0,3.0,1.0e10,1.0e10}; |
---|
675 | double rowUpper[5]={14.0,3.0,3.0,-1.0e10,-1.0e10}; |
---|
676 | double colLower[7]={0.0,0.0,0.0,0.0,0.0,0.0,0.0}; |
---|
677 | double colUpper[7]={100.0,100.0,100.0,100.0,100.0,100.0,100.0}; |
---|
678 | |
---|
679 | // basis 1 |
---|
680 | int rowBasis1[3]={-1,-1,-1}; |
---|
681 | int colBasis1[5]={1,1,-1,-1,1}; |
---|
682 | solution.loadProblem(matrix,colLower,colUpper,objective, |
---|
683 | rowLower,rowUpper); |
---|
684 | int i; |
---|
685 | solution.createStatus(); |
---|
686 | for (i=0;i<3;i++) { |
---|
687 | if (rowBasis1[i]<0) { |
---|
688 | solution.setRowStatus(i,ClpSimplex::atLowerBound); |
---|
689 | } else { |
---|
690 | solution.setRowStatus(i,ClpSimplex::basic); |
---|
691 | } |
---|
692 | } |
---|
693 | for (i=0;i<5;i++) { |
---|
694 | if (colBasis1[i]<0) { |
---|
695 | solution.setColumnStatus(i,ClpSimplex::atLowerBound); |
---|
696 | } else { |
---|
697 | solution.setColumnStatus(i,ClpSimplex::basic); |
---|
698 | } |
---|
699 | } |
---|
700 | solution.setLogLevel(3+4+8+16+32); |
---|
701 | solution.primal(); |
---|
702 | for (i=0;i<3;i++) { |
---|
703 | if (rowBasis1[i]<0) { |
---|
704 | solution.setRowStatus(i,ClpSimplex::atLowerBound); |
---|
705 | } else { |
---|
706 | solution.setRowStatus(i,ClpSimplex::basic); |
---|
707 | } |
---|
708 | } |
---|
709 | for (i=0;i<5;i++) { |
---|
710 | if (colBasis1[i]<0) { |
---|
711 | solution.setColumnStatus(i,ClpSimplex::atLowerBound); |
---|
712 | } else { |
---|
713 | solution.setColumnStatus(i,ClpSimplex::basic); |
---|
714 | } |
---|
715 | } |
---|
716 | // intricate stuff does not work with scaling |
---|
717 | solution.scaling(0); |
---|
718 | #ifndef NDEBUG |
---|
719 | int returnCode = solution.factorize ( ); |
---|
720 | assert(!returnCode); |
---|
721 | #else |
---|
722 | solution.factorize ( ); |
---|
723 | #endif |
---|
724 | const double * colsol = solution.primalColumnSolution(); |
---|
725 | const double * rowsol = solution.primalRowSolution(); |
---|
726 | solution.getSolution(rowsol,colsol); |
---|
727 | #ifndef NDEBUG |
---|
728 | double colsol1[5]={20.0/7.0,3.0,0.0,0.0,23.0/7.0}; |
---|
729 | for (i=0;i<5;i++) { |
---|
730 | assert(eq(colsol[i],colsol1[i])); |
---|
731 | } |
---|
732 | // now feed in again without actually doing factorization |
---|
733 | #endif |
---|
734 | ClpFactorization factorization2 = *solution.factorization(); |
---|
735 | ClpSimplex solution2 = solution; |
---|
736 | solution2.setFactorization(factorization2); |
---|
737 | solution2.createStatus(); |
---|
738 | for (i=0;i<3;i++) { |
---|
739 | if (rowBasis1[i]<0) { |
---|
740 | solution2.setRowStatus(i,ClpSimplex::atLowerBound); |
---|
741 | } else { |
---|
742 | solution2.setRowStatus(i,ClpSimplex::basic); |
---|
743 | } |
---|
744 | } |
---|
745 | for (i=0;i<5;i++) { |
---|
746 | if (colBasis1[i]<0) { |
---|
747 | solution2.setColumnStatus(i,ClpSimplex::atLowerBound); |
---|
748 | } else { |
---|
749 | solution2.setColumnStatus(i,ClpSimplex::basic); |
---|
750 | } |
---|
751 | } |
---|
752 | // intricate stuff does not work with scaling |
---|
753 | solution2.scaling(0); |
---|
754 | solution2.getSolution(rowsol,colsol); |
---|
755 | colsol = solution2.primalColumnSolution(); |
---|
756 | rowsol = solution2.primalRowSolution(); |
---|
757 | for (i=0;i<5;i++) { |
---|
758 | assert(eq(colsol[i],colsol1[i])); |
---|
759 | } |
---|
760 | solution2.setDualBound(0.1); |
---|
761 | solution2.dual(); |
---|
762 | objective[2]=-1.0; |
---|
763 | objective[3]=-0.5; |
---|
764 | objective[4]=10.0; |
---|
765 | solution.dual(); |
---|
766 | for (i=0;i<3;i++) { |
---|
767 | rowLower[i]=-1.0e20; |
---|
768 | colUpper[i+2]=0.0; |
---|
769 | } |
---|
770 | solution.setLogLevel(3); |
---|
771 | solution.dual(); |
---|
772 | double rowObjective[]={1.0,0.5,-10.0}; |
---|
773 | solution.loadProblem(matrix,colLower,colUpper,objective, |
---|
774 | rowLower,rowUpper,rowObjective); |
---|
775 | solution.dual(); |
---|
776 | solution.loadProblem(matrix,colLower,colUpper,objective, |
---|
777 | rowLower,rowUpper,rowObjective); |
---|
778 | solution.primal(); |
---|
779 | } |
---|
780 | #ifndef COIN_NO_CLP_MESSAGE |
---|
781 | { |
---|
782 | CoinMpsIO m; |
---|
783 | std::string fn = dirSample+"exmip1"; |
---|
784 | m.readMps(fn.c_str(),"mps"); |
---|
785 | ClpSimplex solution; |
---|
786 | solution.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(), |
---|
787 | m.getObjCoefficients(), |
---|
788 | m.getRowLower(),m.getRowUpper()); |
---|
789 | solution.dual(); |
---|
790 | // Test event handling |
---|
791 | MyEventHandler handler; |
---|
792 | solution.passInEventHandler(&handler); |
---|
793 | int numberRows=solution.numberRows(); |
---|
794 | // make sure values pass has something to do |
---|
795 | for (int i=0;i<numberRows;i++) |
---|
796 | solution.setRowStatus(i,ClpSimplex::basic); |
---|
797 | solution.primal(1); |
---|
798 | assert (solution.secondaryStatus()==102); // Came out at end of pass |
---|
799 | } |
---|
800 | // Test Message handler |
---|
801 | { |
---|
802 | CoinMpsIO m; |
---|
803 | std::string fn = dirSample+"exmip1"; |
---|
804 | //fn = "Test/subGams4"; |
---|
805 | m.readMps(fn.c_str(),"mps"); |
---|
806 | ClpSimplex model; |
---|
807 | model.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(), |
---|
808 | m.getObjCoefficients(), |
---|
809 | m.getRowLower(),m.getRowUpper()); |
---|
810 | // Message handler |
---|
811 | MyMessageHandler messageHandler(&model); |
---|
812 | std::cout<<"Testing derived message handler"<<std::endl; |
---|
813 | model.passInMessageHandler(&messageHandler); |
---|
814 | model.messagesPointer()->setDetailMessage(1,102); |
---|
815 | model.setFactorizationFrequency(10); |
---|
816 | model.primal(); |
---|
817 | model.primal(0,3); |
---|
818 | model.setObjCoeff(3,-2.9473684210526314); |
---|
819 | model.primal(0,3); |
---|
820 | // Write saved solutions |
---|
821 | int nc = model.getNumCols(); |
---|
822 | int s; |
---|
823 | std::deque<StdVectorDouble> fep = messageHandler.getFeasibleExtremePoints(); |
---|
824 | int numSavedSolutions = fep.size(); |
---|
825 | for ( s=0; s<numSavedSolutions; ++s ) { |
---|
826 | const StdVectorDouble & solnVec = fep[s]; |
---|
827 | for ( int c=0; c<nc; ++c ) { |
---|
828 | if (fabs(solnVec[c])>1.0e-8) |
---|
829 | std::cout <<"Saved Solution: " <<s <<" ColNum: " <<c <<" Value: " <<solnVec[c] <<std::endl; |
---|
830 | } |
---|
831 | } |
---|
832 | // Solve again without scaling |
---|
833 | // and maximize then minimize |
---|
834 | messageHandler.clearFeasibleExtremePoints(); |
---|
835 | model.scaling(0); |
---|
836 | model.setOptimizationDirection(-1); |
---|
837 | model.primal(); |
---|
838 | model.setOptimizationDirection(1); |
---|
839 | model.primal(); |
---|
840 | fep = messageHandler.getFeasibleExtremePoints(); |
---|
841 | numSavedSolutions = fep.size(); |
---|
842 | for ( s=0; s<numSavedSolutions; ++s ) { |
---|
843 | const StdVectorDouble & solnVec = fep[s]; |
---|
844 | for ( int c=0; c<nc; ++c ) { |
---|
845 | if (fabs(solnVec[c])>1.0e-8) |
---|
846 | std::cout <<"Saved Solution: " <<s <<" ColNum: " <<c <<" Value: " <<solnVec[c] <<std::endl; |
---|
847 | } |
---|
848 | } |
---|
849 | } |
---|
850 | #endif |
---|
851 | // Test dual ranging |
---|
852 | { |
---|
853 | CoinMpsIO m; |
---|
854 | std::string fn = dirSample+"exmip1"; |
---|
855 | m.readMps(fn.c_str(),"mps"); |
---|
856 | ClpSimplex model; |
---|
857 | model.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(), |
---|
858 | m.getObjCoefficients(), |
---|
859 | m.getRowLower(),m.getRowUpper()); |
---|
860 | model.primal(); |
---|
861 | int which[13] = {0,1,2,3,4,5,6,7,8,9,10,11,12}; |
---|
862 | double costIncrease[13]; |
---|
863 | int sequenceIncrease[13]; |
---|
864 | double costDecrease[13]; |
---|
865 | int sequenceDecrease[13]; |
---|
866 | // ranging |
---|
867 | model.dualRanging(13,which,costIncrease,sequenceIncrease, |
---|
868 | costDecrease,sequenceDecrease); |
---|
869 | int i; |
---|
870 | for ( i=0;i<13;i++) |
---|
871 | printf("%d increase %g %d, decrease %g %d\n", |
---|
872 | i,costIncrease[i],sequenceIncrease[i], |
---|
873 | costDecrease[i],sequenceDecrease[i]); |
---|
874 | assert (fabs(costDecrease[3])<1.0e-4); |
---|
875 | assert (fabs(costIncrease[7]-1.0)<1.0e-4); |
---|
876 | model.setOptimizationDirection(-1); |
---|
877 | { |
---|
878 | int j; |
---|
879 | double * obj = model.objective(); |
---|
880 | int n=model.numberColumns(); |
---|
881 | for (j=0;j<n;j++) |
---|
882 | obj[j] *= -1.0; |
---|
883 | } |
---|
884 | double costIncrease2[13]; |
---|
885 | int sequenceIncrease2[13]; |
---|
886 | double costDecrease2[13]; |
---|
887 | int sequenceDecrease2[13]; |
---|
888 | // ranging |
---|
889 | model.dualRanging(13,which,costIncrease2,sequenceIncrease2, |
---|
890 | costDecrease2,sequenceDecrease2); |
---|
891 | for (i=0;i<13;i++) { |
---|
892 | assert (fabs(costIncrease[i]-costDecrease2[i])<1.0e-6); |
---|
893 | assert (fabs(costDecrease[i]-costIncrease2[i])<1.0e-6); |
---|
894 | assert (sequenceIncrease[i]==sequenceDecrease2[i]); |
---|
895 | assert (sequenceDecrease[i]==sequenceIncrease2[i]); |
---|
896 | } |
---|
897 | // Now delete all rows and see what happens |
---|
898 | model.deleteRows(model.numberRows(),which); |
---|
899 | model.primal(); |
---|
900 | // ranging |
---|
901 | if (!model.dualRanging(8,which,costIncrease,sequenceIncrease, |
---|
902 | costDecrease,sequenceDecrease)) { |
---|
903 | for (i=0;i<8;i++) { |
---|
904 | printf("%d increase %g %d, decrease %g %d\n", |
---|
905 | i,costIncrease[i],sequenceIncrease[i], |
---|
906 | costDecrease[i],sequenceDecrease[i]); |
---|
907 | } |
---|
908 | } |
---|
909 | } |
---|
910 | // Test primal ranging |
---|
911 | { |
---|
912 | CoinMpsIO m; |
---|
913 | std::string fn = dirSample+"exmip1"; |
---|
914 | m.readMps(fn.c_str(),"mps"); |
---|
915 | ClpSimplex model; |
---|
916 | model.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(), |
---|
917 | m.getObjCoefficients(), |
---|
918 | m.getRowLower(),m.getRowUpper()); |
---|
919 | model.primal(); |
---|
920 | int which[13] = {0,1,2,3,4,5,6,7,8,9,10,11,12}; |
---|
921 | double valueIncrease[13]; |
---|
922 | int sequenceIncrease[13]; |
---|
923 | double valueDecrease[13]; |
---|
924 | int sequenceDecrease[13]; |
---|
925 | // ranging |
---|
926 | model.primalRanging(13,which,valueIncrease,sequenceIncrease, |
---|
927 | valueDecrease,sequenceDecrease); |
---|
928 | int i; |
---|
929 | for ( i=0;i<13;i++) |
---|
930 | printf("%d increase %g %d, decrease %g %d\n", |
---|
931 | i,valueIncrease[i],sequenceIncrease[i], |
---|
932 | valueDecrease[i],sequenceDecrease[i]); |
---|
933 | assert (fabs(valueIncrease[3]-0.642857)<1.0e-4); |
---|
934 | assert (fabs(valueIncrease[8]-2.95113)<1.0e-4); |
---|
935 | #if 0 |
---|
936 | // out until I find optimization bug |
---|
937 | // Test parametrics |
---|
938 | ClpSimplexOther * model2 = (ClpSimplexOther *) (&model); |
---|
939 | double rhs[]={ 1.0,2.0,3.0,4.0,5.0}; |
---|
940 | double endingTheta=1.0; |
---|
941 | model2->scaling(0); |
---|
942 | model2->setLogLevel(63); |
---|
943 | model2->parametrics(0.0,endingTheta,0.1, |
---|
944 | NULL,NULL,rhs,rhs,NULL); |
---|
945 | #endif |
---|
946 | } |
---|
947 | // Test binv etc |
---|
948 | { |
---|
949 | /* |
---|
950 | Wolsey : Page 130 |
---|
951 | max 4x1 - x2 |
---|
952 | 7x1 - 2x2 <= 14 |
---|
953 | x2 <= 3 |
---|
954 | 2x1 - 2x2 <= 3 |
---|
955 | x1 in Z+, x2 >= 0 |
---|
956 | |
---|
957 | note slacks are -1 in Clp so signs may be different |
---|
958 | */ |
---|
959 | |
---|
960 | int n_cols = 2; |
---|
961 | int n_rows = 3; |
---|
962 | |
---|
963 | double obj[2] = {-4.0, 1.0}; |
---|
964 | double collb[2] = {0.0, 0.0}; |
---|
965 | double colub[2] = {COIN_DBL_MAX, COIN_DBL_MAX}; |
---|
966 | double rowlb[3] = {-COIN_DBL_MAX, -COIN_DBL_MAX, -COIN_DBL_MAX}; |
---|
967 | double rowub[3] = {14.0, 3.0, 3.0}; |
---|
968 | |
---|
969 | int rowIndices[5] = {0, 2, 0, 1, 2}; |
---|
970 | int colIndices[5] = {0, 0, 1, 1, 1}; |
---|
971 | double elements[5] = {7.0, 2.0, -2.0, 1.0, -2.0}; |
---|
972 | CoinPackedMatrix M(true, rowIndices, colIndices, elements, 5); |
---|
973 | |
---|
974 | ClpSimplex model; |
---|
975 | model.loadProblem(M, collb, colub, obj, rowlb, rowub); |
---|
976 | model.dual(0,1); // keep factorization |
---|
977 | |
---|
978 | //check that the tableau matches wolsey (B-1 A) |
---|
979 | // slacks in second part of binvA |
---|
980 | double * binvA = reinterpret_cast<double*> (malloc((n_cols+n_rows) * sizeof(double))); |
---|
981 | |
---|
982 | printf("B-1 A by row\n"); |
---|
983 | int i; |
---|
984 | for( i = 0; i < n_rows; i++){ |
---|
985 | model.getBInvARow(i, binvA,binvA+n_cols); |
---|
986 | printf("row: %d -> ",i); |
---|
987 | for(int j=0; j < n_cols+n_rows; j++){ |
---|
988 | printf("%g, ", binvA[j]); |
---|
989 | } |
---|
990 | printf("\n"); |
---|
991 | } |
---|
992 | // See if can re-use factorization AND arrays |
---|
993 | model.primal(0,3+4); // keep factorization |
---|
994 | // And do by column |
---|
995 | printf("B-1 A by column\n"); |
---|
996 | for( i = 0; i < n_rows+n_cols; i++){ |
---|
997 | model.getBInvACol(i, binvA); |
---|
998 | printf("column: %d -> ",i); |
---|
999 | for(int j=0; j < n_rows; j++){ |
---|
1000 | printf("%g, ", binvA[j]); |
---|
1001 | } |
---|
1002 | printf("\n"); |
---|
1003 | } |
---|
1004 | /* Do twice - |
---|
1005 | without and with scaling |
---|
1006 | */ |
---|
1007 | // set scaling off |
---|
1008 | model.scaling(0); |
---|
1009 | for (int iPass=0;iPass<2;iPass++) { |
---|
1010 | model.primal(0,3+4); // keep factorization |
---|
1011 | const double * rowScale = model.rowScale(); |
---|
1012 | const double * columnScale = model.columnScale(); |
---|
1013 | if (!iPass) |
---|
1014 | assert (!rowScale); |
---|
1015 | else |
---|
1016 | assert (rowScale); // only true for this example |
---|
1017 | /* has to be exactly correct as in OsiClpsolverInterface.cpp |
---|
1018 | (also redo each pass as may change |
---|
1019 | */ |
---|
1020 | printf("B-1 A"); |
---|
1021 | for( i = 0; i < n_rows; i++){ |
---|
1022 | model.getBInvARow(i, binvA,binvA+n_cols); |
---|
1023 | printf("\nrow: %d -> ",i); |
---|
1024 | int j; |
---|
1025 | // First columns |
---|
1026 | for(j=0; j < n_cols; j++){ |
---|
1027 | if (binvA[j]) { |
---|
1028 | printf("(%d %g), ", j, binvA[j]); |
---|
1029 | } |
---|
1030 | } |
---|
1031 | // now rows |
---|
1032 | for(j=0; j < n_rows; j++){ |
---|
1033 | if (binvA[j+n_cols]) { |
---|
1034 | printf("(%d %g), ", j+n_cols, binvA[j+n_cols]); |
---|
1035 | } |
---|
1036 | } |
---|
1037 | } |
---|
1038 | printf("\n"); |
---|
1039 | printf("And by column (trickier)"); |
---|
1040 | const int * pivotVariable = model.pivotVariable(); |
---|
1041 | for( i = 0; i < n_cols+n_rows; i++){ |
---|
1042 | model.getBInvACol(i, binvA); |
---|
1043 | printf("\ncolumn: %d -> ",i); |
---|
1044 | for(int j=0; j < n_rows; j++){ |
---|
1045 | if (binvA[j]) { |
---|
1046 | // need to know pivot variable for +1/-1 (slack) and row/column scaling |
---|
1047 | int pivot = pivotVariable[j]; |
---|
1048 | if (pivot<n_cols) { |
---|
1049 | // scaled coding is in just in case |
---|
1050 | if (!columnScale) { |
---|
1051 | printf("(%d %g), ", j, binvA[j]); |
---|
1052 | } else { |
---|
1053 | printf("(%d %g), ", j, binvA[j]*columnScale[pivot]); |
---|
1054 | } |
---|
1055 | } else { |
---|
1056 | if (!rowScale) { |
---|
1057 | printf("(%d %g), ", j, binvA[j]); |
---|
1058 | } else { |
---|
1059 | printf("(%d %g), ", j, binvA[j]/rowScale[pivot-n_cols]); |
---|
1060 | } |
---|
1061 | } |
---|
1062 | } |
---|
1063 | } |
---|
1064 | } |
---|
1065 | printf("\n"); |
---|
1066 | printf("binvrow"); |
---|
1067 | for( i = 0; i < n_rows; i++){ |
---|
1068 | model.getBInvRow(i, binvA); |
---|
1069 | printf("\nrow: %d -> ",i); |
---|
1070 | int j; |
---|
1071 | for (j=0;j<n_rows;j++) { |
---|
1072 | if (binvA[j]) |
---|
1073 | printf("(%d %g), ", j, binvA[j]); |
---|
1074 | } |
---|
1075 | } |
---|
1076 | printf("\n"); |
---|
1077 | printf("And by column "); |
---|
1078 | for( i = 0; i < n_rows; i++){ |
---|
1079 | model.getBInvCol(i, binvA); |
---|
1080 | printf("\ncol: %d -> ",i); |
---|
1081 | int j; |
---|
1082 | for (j=0;j<n_rows;j++) { |
---|
1083 | if (binvA[j]) |
---|
1084 | printf("(%d %g), ", j, binvA[j]); |
---|
1085 | } |
---|
1086 | } |
---|
1087 | printf("\n"); |
---|
1088 | // now deal with next pass |
---|
1089 | if (!iPass) { |
---|
1090 | // get scaling for testing |
---|
1091 | model.scaling(1); |
---|
1092 | } |
---|
1093 | } |
---|
1094 | free(binvA); |
---|
1095 | model.setColUpper(1,2.0); |
---|
1096 | model.dual(0,2+4); // use factorization and arrays |
---|
1097 | model.dual(0,2); // hopefully will not use factorization |
---|
1098 | model.primal(0,3+4); // keep factorization |
---|
1099 | // but say basis has changed |
---|
1100 | model.setWhatsChanged(model.whatsChanged()&(~512)); |
---|
1101 | model.dual(0,2); // hopefully will not use factorization |
---|
1102 | } |
---|
1103 | // test steepest edge |
---|
1104 | { |
---|
1105 | CoinMpsIO m; |
---|
1106 | std::string fn = dirSample+"finnis"; |
---|
1107 | int returnCode = m.readMps(fn.c_str(),"mps"); |
---|
1108 | if (returnCode) { |
---|
1109 | // probable cause is that gz not there |
---|
1110 | fprintf(stderr,"Unable to open finnis.mps in %s\n", dirSample.c_str()); |
---|
1111 | fprintf(stderr,"Most probable cause is finnis.mps is gzipped i.e. finnis.mps.gz and libz has not been activated\n"); |
---|
1112 | fprintf(stderr,"Either gunzip files or edit Makefiles/Makefile.location to get libz\n"); |
---|
1113 | exit(999); |
---|
1114 | } |
---|
1115 | ClpModel model; |
---|
1116 | model.loadProblem(*m.getMatrixByCol(),m.getColLower(), |
---|
1117 | m.getColUpper(), |
---|
1118 | m.getObjCoefficients(), |
---|
1119 | m.getRowLower(),m.getRowUpper()); |
---|
1120 | ClpSimplex solution(model); |
---|
1121 | |
---|
1122 | solution.scaling(1); |
---|
1123 | solution.setDualBound(1.0e8); |
---|
1124 | //solution.factorization()->maximumPivots(1); |
---|
1125 | //solution.setLogLevel(3); |
---|
1126 | solution.setDualTolerance(1.0e-7); |
---|
1127 | // set objective sense, |
---|
1128 | ClpDualRowSteepest steep; |
---|
1129 | solution.setDualRowPivotAlgorithm(steep); |
---|
1130 | solution.setDblParam(ClpObjOffset,m.objectiveOffset()); |
---|
1131 | solution.dual(); |
---|
1132 | } |
---|
1133 | // test normal solution |
---|
1134 | { |
---|
1135 | CoinMpsIO m; |
---|
1136 | std::string fn = dirSample+"afiro"; |
---|
1137 | m.readMps(fn.c_str(),"mps"); |
---|
1138 | ClpSimplex solution; |
---|
1139 | ClpModel model; |
---|
1140 | // do twice - without and with scaling |
---|
1141 | int iPass; |
---|
1142 | for (iPass=0;iPass<2;iPass++) { |
---|
1143 | // explicit row objective for testing |
---|
1144 | int nr = m.getNumRows(); |
---|
1145 | double * rowObj = new double[nr]; |
---|
1146 | CoinFillN(rowObj,nr,0.0); |
---|
1147 | model.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(), |
---|
1148 | m.getObjCoefficients(), |
---|
1149 | m.getRowLower(),m.getRowUpper(),rowObj); |
---|
1150 | delete [] rowObj; |
---|
1151 | solution = ClpSimplex(model); |
---|
1152 | if (iPass) { |
---|
1153 | solution.scaling(); |
---|
1154 | } |
---|
1155 | solution.dual(); |
---|
1156 | solution.dual(); |
---|
1157 | // test optimal |
---|
1158 | assert (solution.status()==0); |
---|
1159 | int numberColumns = solution.numberColumns(); |
---|
1160 | int numberRows = solution.numberRows(); |
---|
1161 | CoinPackedVector colsol(numberColumns,solution.primalColumnSolution()); |
---|
1162 | double * objective = solution.objective(); |
---|
1163 | #ifndef NDEBUG |
---|
1164 | double objValue = colsol.dotProduct(objective); |
---|
1165 | #endif |
---|
1166 | CoinRelFltEq eq(1.0e-8); |
---|
1167 | assert(eq(objValue,-4.6475314286e+02)); |
---|
1168 | solution.dual(); |
---|
1169 | assert(eq(solution.objectiveValue(),-4.6475314286e+02)); |
---|
1170 | double * lower = solution.columnLower(); |
---|
1171 | double * upper = solution.columnUpper(); |
---|
1172 | double * sol = solution.primalColumnSolution(); |
---|
1173 | double * result = new double[numberColumns]; |
---|
1174 | CoinFillN ( result, numberColumns,0.0); |
---|
1175 | solution.matrix()->transposeTimes(solution.dualRowSolution(), result); |
---|
1176 | int iRow , iColumn; |
---|
1177 | // see if feasible and dual feasible |
---|
1178 | for (iColumn=0;iColumn<numberColumns;iColumn++) { |
---|
1179 | double value = sol[iColumn]; |
---|
1180 | assert(value<upper[iColumn]+1.0e-8); |
---|
1181 | assert(value>lower[iColumn]-1.0e-8); |
---|
1182 | value = objective[iColumn]-result[iColumn]; |
---|
1183 | assert (value>-1.0e-5); |
---|
1184 | if (sol[iColumn]>1.0e-5) |
---|
1185 | assert (value<1.0e-5); |
---|
1186 | } |
---|
1187 | delete [] result; |
---|
1188 | result = new double[numberRows]; |
---|
1189 | CoinFillN ( result, numberRows,0.0); |
---|
1190 | solution.matrix()->times(colsol, result); |
---|
1191 | lower = solution.rowLower(); |
---|
1192 | upper = solution.rowUpper(); |
---|
1193 | sol = solution.primalRowSolution(); |
---|
1194 | #ifndef NDEBUG |
---|
1195 | for (iRow=0;iRow<numberRows;iRow++) { |
---|
1196 | double value = result[iRow]; |
---|
1197 | assert(eq(value,sol[iRow])); |
---|
1198 | assert(value<upper[iRow]+1.0e-8); |
---|
1199 | assert(value>lower[iRow]-1.0e-8); |
---|
1200 | } |
---|
1201 | #endif |
---|
1202 | delete [] result; |
---|
1203 | // test row objective |
---|
1204 | double * rowObjective = solution.rowObjective(); |
---|
1205 | CoinDisjointCopyN(solution.dualRowSolution(),numberRows,rowObjective); |
---|
1206 | CoinDisjointCopyN(solution.dualColumnSolution(),numberColumns,objective); |
---|
1207 | // this sets up all slack basis |
---|
1208 | solution.createStatus(); |
---|
1209 | solution.dual(); |
---|
1210 | CoinFillN(rowObjective,numberRows,0.0); |
---|
1211 | CoinDisjointCopyN(m.getObjCoefficients(),numberColumns,objective); |
---|
1212 | solution.dual(); |
---|
1213 | } |
---|
1214 | } |
---|
1215 | // test unbounded |
---|
1216 | { |
---|
1217 | CoinMpsIO m; |
---|
1218 | std::string fn = dirSample+"brandy"; |
---|
1219 | m.readMps(fn.c_str(),"mps"); |
---|
1220 | ClpSimplex solution; |
---|
1221 | // do twice - without and with scaling |
---|
1222 | int iPass; |
---|
1223 | for (iPass=0;iPass<2;iPass++) { |
---|
1224 | solution.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(), |
---|
1225 | m.getObjCoefficients(), |
---|
1226 | m.getRowLower(),m.getRowUpper()); |
---|
1227 | if (iPass) |
---|
1228 | solution.scaling(); |
---|
1229 | solution.setOptimizationDirection(-1); |
---|
1230 | // test unbounded and ray |
---|
1231 | #ifdef DUAL |
---|
1232 | solution.setDualBound(100.0); |
---|
1233 | solution.dual(); |
---|
1234 | #else |
---|
1235 | solution.primal(); |
---|
1236 | #endif |
---|
1237 | assert (solution.status()==2); |
---|
1238 | int numberColumns = solution.numberColumns(); |
---|
1239 | int numberRows = solution.numberRows(); |
---|
1240 | double * lower = solution.columnLower(); |
---|
1241 | double * upper = solution.columnUpper(); |
---|
1242 | double * sol = solution.primalColumnSolution(); |
---|
1243 | double * ray = solution.unboundedRay(); |
---|
1244 | double * objective = solution.objective(); |
---|
1245 | double objChange=0.0; |
---|
1246 | int iRow , iColumn; |
---|
1247 | // make sure feasible and columns form ray |
---|
1248 | for (iColumn=0;iColumn<numberColumns;iColumn++) { |
---|
1249 | double value = sol[iColumn]; |
---|
1250 | assert(value<upper[iColumn]+1.0e-8); |
---|
1251 | assert(value>lower[iColumn]-1.0e-8); |
---|
1252 | value = ray[iColumn]; |
---|
1253 | if (value>0.0) |
---|
1254 | assert(upper[iColumn]>1.0e30); |
---|
1255 | else if (value<0.0) |
---|
1256 | assert(lower[iColumn]<-1.0e30); |
---|
1257 | objChange += value*objective[iColumn]; |
---|
1258 | } |
---|
1259 | // make sure increasing objective |
---|
1260 | assert(objChange>0.0); |
---|
1261 | double * result = new double[numberRows]; |
---|
1262 | CoinFillN ( result, numberRows,0.0); |
---|
1263 | solution.matrix()->times(sol, result); |
---|
1264 | lower = solution.rowLower(); |
---|
1265 | upper = solution.rowUpper(); |
---|
1266 | sol = solution.primalRowSolution(); |
---|
1267 | #ifndef NDEBUG |
---|
1268 | for (iRow=0;iRow<numberRows;iRow++) { |
---|
1269 | double value = result[iRow]; |
---|
1270 | assert(eq(value,sol[iRow])); |
---|
1271 | assert(value<upper[iRow]+2.0e-8); |
---|
1272 | assert(value>lower[iRow]-2.0e-8); |
---|
1273 | } |
---|
1274 | #endif |
---|
1275 | CoinFillN ( result, numberRows,0.0); |
---|
1276 | solution.matrix()->times(ray, result); |
---|
1277 | // there may be small differences (especially if scaled) |
---|
1278 | for (iRow=0;iRow<numberRows;iRow++) { |
---|
1279 | double value = result[iRow]; |
---|
1280 | if (value>1.0e-8) |
---|
1281 | assert(upper[iRow]>1.0e30); |
---|
1282 | else if (value<-1.0e-8) |
---|
1283 | assert(lower[iRow]<-1.0e30); |
---|
1284 | } |
---|
1285 | delete [] result; |
---|
1286 | delete [] ray; |
---|
1287 | } |
---|
1288 | } |
---|
1289 | // test infeasible |
---|
1290 | { |
---|
1291 | CoinMpsIO m; |
---|
1292 | std::string fn = dirSample+"brandy"; |
---|
1293 | m.readMps(fn.c_str(),"mps"); |
---|
1294 | ClpSimplex solution; |
---|
1295 | // do twice - without and with scaling |
---|
1296 | int iPass; |
---|
1297 | for (iPass=0;iPass<2;iPass++) { |
---|
1298 | solution.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(), |
---|
1299 | m.getObjCoefficients(), |
---|
1300 | m.getRowLower(),m.getRowUpper()); |
---|
1301 | if (iPass) |
---|
1302 | solution.scaling(); |
---|
1303 | // test infeasible and ray |
---|
1304 | solution.columnUpper()[0]=0.0; |
---|
1305 | #ifdef DUAL |
---|
1306 | solution.setDualBound(100.0); |
---|
1307 | solution.dual(); |
---|
1308 | #else |
---|
1309 | solution.primal(); |
---|
1310 | #endif |
---|
1311 | assert (solution.status()==1); |
---|
1312 | int numberColumns = solution.numberColumns(); |
---|
1313 | int numberRows = solution.numberRows(); |
---|
1314 | double * lower = solution.rowLower(); |
---|
1315 | double * upper = solution.rowUpper(); |
---|
1316 | double * ray = solution.infeasibilityRay(); |
---|
1317 | assert(ray); |
---|
1318 | // construct proof of infeasibility |
---|
1319 | int iRow , iColumn; |
---|
1320 | double lo=0.0,up=0.0; |
---|
1321 | int nl=0,nu=0; |
---|
1322 | for (iRow=0;iRow<numberRows;iRow++) { |
---|
1323 | if (lower[iRow]>-1.0e20) { |
---|
1324 | lo += ray[iRow]*lower[iRow]; |
---|
1325 | } else { |
---|
1326 | if (ray[iRow]>1.0e-8) |
---|
1327 | nl++; |
---|
1328 | } |
---|
1329 | if (upper[iRow]<1.0e20) { |
---|
1330 | up += ray[iRow]*upper[iRow]; |
---|
1331 | } else { |
---|
1332 | if (ray[iRow]>1.0e-8) |
---|
1333 | nu++; |
---|
1334 | } |
---|
1335 | } |
---|
1336 | if (nl) |
---|
1337 | lo=-1.0e100; |
---|
1338 | if (nu) |
---|
1339 | up=1.0e100; |
---|
1340 | double * result = new double[numberColumns]; |
---|
1341 | double lo2=0.0,up2=0.0; |
---|
1342 | CoinFillN ( result, numberColumns,0.0); |
---|
1343 | solution.matrix()->transposeTimes(ray, result); |
---|
1344 | lower = solution.columnLower(); |
---|
1345 | upper = solution.columnUpper(); |
---|
1346 | nl=nu=0; |
---|
1347 | for (iColumn=0;iColumn<numberColumns;iColumn++) { |
---|
1348 | if (result[iColumn]>1.0e-8) { |
---|
1349 | if (lower[iColumn]>-1.0e20) |
---|
1350 | lo2 += result[iColumn]*lower[iColumn]; |
---|
1351 | else |
---|
1352 | nl++; |
---|
1353 | if (upper[iColumn]<1.0e20) |
---|
1354 | up2 += result[iColumn]*upper[iColumn]; |
---|
1355 | else |
---|
1356 | nu++; |
---|
1357 | } else if (result[iColumn]<-1.0e-8) { |
---|
1358 | if (lower[iColumn]>-1.0e20) |
---|
1359 | up2 += result[iColumn]*lower[iColumn]; |
---|
1360 | else |
---|
1361 | nu++; |
---|
1362 | if (upper[iColumn]<1.0e20) |
---|
1363 | lo2 += result[iColumn]*upper[iColumn]; |
---|
1364 | else |
---|
1365 | nl++; |
---|
1366 | } |
---|
1367 | } |
---|
1368 | if (nl) |
---|
1369 | lo2=-1.0e100; |
---|
1370 | if (nu) |
---|
1371 | up2=1.0e100; |
---|
1372 | // make sure inconsistency |
---|
1373 | assert(lo2>up||up2<lo); |
---|
1374 | delete [] result; |
---|
1375 | delete [] ray; |
---|
1376 | } |
---|
1377 | } |
---|
1378 | // test delete and add |
---|
1379 | { |
---|
1380 | CoinMpsIO m; |
---|
1381 | std::string fn = dirSample+"brandy"; |
---|
1382 | m.readMps(fn.c_str(),"mps"); |
---|
1383 | ClpSimplex solution; |
---|
1384 | solution.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(), |
---|
1385 | m.getObjCoefficients(), |
---|
1386 | m.getRowLower(),m.getRowUpper()); |
---|
1387 | solution.dual(); |
---|
1388 | CoinRelFltEq eq(1.0e-8); |
---|
1389 | assert(eq(solution.objectiveValue(),1.5185098965e+03)); |
---|
1390 | |
---|
1391 | int numberColumns = solution.numberColumns(); |
---|
1392 | int numberRows = solution.numberRows(); |
---|
1393 | double * saveObj = new double [numberColumns]; |
---|
1394 | double * saveLower = new double[numberRows+numberColumns]; |
---|
1395 | double * saveUpper = new double[numberRows+numberColumns]; |
---|
1396 | int * which = new int [numberRows+numberColumns]; |
---|
1397 | |
---|
1398 | int numberElements = m.getMatrixByCol()->getNumElements(); |
---|
1399 | int * starts = new int[numberRows+numberColumns]; |
---|
1400 | int * index = new int[numberElements]; |
---|
1401 | double * element = new double[numberElements]; |
---|
1402 | |
---|
1403 | const CoinBigIndex * startM; |
---|
1404 | const int * lengthM; |
---|
1405 | const int * indexM; |
---|
1406 | const double * elementM; |
---|
1407 | |
---|
1408 | int n,nel; |
---|
1409 | |
---|
1410 | // delete non basic columns |
---|
1411 | n=0; |
---|
1412 | nel=0; |
---|
1413 | int iRow , iColumn; |
---|
1414 | const double * lower = m.getColLower(); |
---|
1415 | const double * upper = m.getColUpper(); |
---|
1416 | const double * objective = m.getObjCoefficients(); |
---|
1417 | startM = m.getMatrixByCol()->getVectorStarts(); |
---|
1418 | lengthM = m.getMatrixByCol()->getVectorLengths(); |
---|
1419 | indexM = m.getMatrixByCol()->getIndices(); |
---|
1420 | elementM = m.getMatrixByCol()->getElements(); |
---|
1421 | starts[0]=0; |
---|
1422 | for (iColumn=0;iColumn<numberColumns;iColumn++) { |
---|
1423 | if (solution.getColumnStatus(iColumn)!=ClpSimplex::basic) { |
---|
1424 | saveObj[n]=objective[iColumn]; |
---|
1425 | saveLower[n]=lower[iColumn]; |
---|
1426 | saveUpper[n]=upper[iColumn]; |
---|
1427 | int j; |
---|
1428 | for (j=startM[iColumn];j<startM[iColumn]+lengthM[iColumn];j++) { |
---|
1429 | index[nel]=indexM[j]; |
---|
1430 | element[nel++]=elementM[j]; |
---|
1431 | } |
---|
1432 | which[n++]=iColumn; |
---|
1433 | starts[n]=nel; |
---|
1434 | } |
---|
1435 | } |
---|
1436 | solution.deleteColumns(n,which); |
---|
1437 | solution.dual(); |
---|
1438 | // Put back |
---|
1439 | solution.addColumns(n,saveLower,saveUpper,saveObj, |
---|
1440 | starts,index,element); |
---|
1441 | solution.dual(); |
---|
1442 | assert(eq(solution.objectiveValue(),1.5185098965e+03)); |
---|
1443 | // Delete all columns and add back |
---|
1444 | n=0; |
---|
1445 | nel=0; |
---|
1446 | starts[0]=0; |
---|
1447 | lower = m.getColLower(); |
---|
1448 | upper = m.getColUpper(); |
---|
1449 | objective = m.getObjCoefficients(); |
---|
1450 | for (iColumn=0;iColumn<numberColumns;iColumn++) { |
---|
1451 | saveObj[n]=objective[iColumn]; |
---|
1452 | saveLower[n]=lower[iColumn]; |
---|
1453 | saveUpper[n]=upper[iColumn]; |
---|
1454 | int j; |
---|
1455 | for (j=startM[iColumn];j<startM[iColumn]+lengthM[iColumn];j++) { |
---|
1456 | index[nel]=indexM[j]; |
---|
1457 | element[nel++]=elementM[j]; |
---|
1458 | } |
---|
1459 | which[n++]=iColumn; |
---|
1460 | starts[n]=nel; |
---|
1461 | } |
---|
1462 | solution.deleteColumns(n,which); |
---|
1463 | solution.dual(); |
---|
1464 | // Put back |
---|
1465 | solution.addColumns(n,saveLower,saveUpper,saveObj, |
---|
1466 | starts,index,element); |
---|
1467 | solution.dual(); |
---|
1468 | assert(eq(solution.objectiveValue(),1.5185098965e+03)); |
---|
1469 | |
---|
1470 | // reload with original |
---|
1471 | solution.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(), |
---|
1472 | m.getObjCoefficients(), |
---|
1473 | m.getRowLower(),m.getRowUpper()); |
---|
1474 | // delete half rows |
---|
1475 | n=0; |
---|
1476 | nel=0; |
---|
1477 | lower = m.getRowLower(); |
---|
1478 | upper = m.getRowUpper(); |
---|
1479 | startM = m.getMatrixByRow()->getVectorStarts(); |
---|
1480 | lengthM = m.getMatrixByRow()->getVectorLengths(); |
---|
1481 | indexM = m.getMatrixByRow()->getIndices(); |
---|
1482 | elementM = m.getMatrixByRow()->getElements(); |
---|
1483 | starts[0]=0; |
---|
1484 | for (iRow=0;iRow<numberRows;iRow++) { |
---|
1485 | if ((iRow&1)==0) { |
---|
1486 | saveLower[n]=lower[iRow]; |
---|
1487 | saveUpper[n]=upper[iRow]; |
---|
1488 | int j; |
---|
1489 | for (j=startM[iRow];j<startM[iRow]+lengthM[iRow];j++) { |
---|
1490 | index[nel]=indexM[j]; |
---|
1491 | element[nel++]=elementM[j]; |
---|
1492 | } |
---|
1493 | which[n++]=iRow; |
---|
1494 | starts[n]=nel; |
---|
1495 | } |
---|
1496 | } |
---|
1497 | solution.deleteRows(n,which); |
---|
1498 | solution.dual(); |
---|
1499 | // Put back |
---|
1500 | solution.addRows(n,saveLower,saveUpper, |
---|
1501 | starts,index,element); |
---|
1502 | solution.dual(); |
---|
1503 | assert(eq(solution.objectiveValue(),1.5185098965e+03)); |
---|
1504 | solution.writeMps("yy.mps"); |
---|
1505 | // Delete all rows |
---|
1506 | n=0; |
---|
1507 | nel=0; |
---|
1508 | lower = m.getRowLower(); |
---|
1509 | upper = m.getRowUpper(); |
---|
1510 | starts[0]=0; |
---|
1511 | for (iRow=0;iRow<numberRows;iRow++) { |
---|
1512 | saveLower[n]=lower[iRow]; |
---|
1513 | saveUpper[n]=upper[iRow]; |
---|
1514 | int j; |
---|
1515 | for (j=startM[iRow];j<startM[iRow]+lengthM[iRow];j++) { |
---|
1516 | index[nel]=indexM[j]; |
---|
1517 | element[nel++]=elementM[j]; |
---|
1518 | } |
---|
1519 | which[n++]=iRow; |
---|
1520 | starts[n]=nel; |
---|
1521 | } |
---|
1522 | solution.deleteRows(n,which); |
---|
1523 | solution.dual(); |
---|
1524 | // Put back |
---|
1525 | solution.addRows(n,saveLower,saveUpper, |
---|
1526 | starts,index,element); |
---|
1527 | solution.dual(); |
---|
1528 | solution.writeMps("xx.mps"); |
---|
1529 | assert(eq(solution.objectiveValue(),1.5185098965e+03)); |
---|
1530 | // Zero out status array to give some interest |
---|
1531 | memset(solution.statusArray()+numberColumns,0,numberRows); |
---|
1532 | solution.primal(1); |
---|
1533 | assert(eq(solution.objectiveValue(),1.5185098965e+03)); |
---|
1534 | // Delete all columns and rows |
---|
1535 | n=0; |
---|
1536 | for (iColumn=0;iColumn<numberColumns;iColumn++) { |
---|
1537 | which[n++]=iColumn; |
---|
1538 | starts[n]=nel; |
---|
1539 | } |
---|
1540 | solution.deleteColumns(n,which); |
---|
1541 | n=0; |
---|
1542 | for (iRow=0;iRow<numberRows;iRow++) { |
---|
1543 | which[n++]=iRow; |
---|
1544 | starts[n]=nel; |
---|
1545 | } |
---|
1546 | solution.deleteRows(n,which); |
---|
1547 | |
---|
1548 | delete [] saveObj; |
---|
1549 | delete [] saveLower; |
---|
1550 | delete [] saveUpper; |
---|
1551 | delete [] which; |
---|
1552 | delete [] starts; |
---|
1553 | delete [] index; |
---|
1554 | delete [] element; |
---|
1555 | } |
---|
1556 | #if 1 |
---|
1557 | // Test barrier |
---|
1558 | { |
---|
1559 | CoinMpsIO m; |
---|
1560 | std::string fn = dirSample+"exmip1"; |
---|
1561 | m.readMps(fn.c_str(),"mps"); |
---|
1562 | ClpInterior solution; |
---|
1563 | solution.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(), |
---|
1564 | m.getObjCoefficients(), |
---|
1565 | m.getRowLower(),m.getRowUpper()); |
---|
1566 | solution.primalDual(); |
---|
1567 | } |
---|
1568 | #endif |
---|
1569 | // test network |
---|
1570 | #define QUADRATIC |
---|
1571 | if (1) { |
---|
1572 | std::string fn = dirSample+"input.130"; |
---|
1573 | int numberColumns; |
---|
1574 | int numberRows; |
---|
1575 | |
---|
1576 | FILE * fp = fopen(fn.c_str(),"r"); |
---|
1577 | if (!fp) { |
---|
1578 | // Try in Data/Sample |
---|
1579 | fn = "Data/Sample/input.130"; |
---|
1580 | fp = fopen(fn.c_str(),"r"); |
---|
1581 | } |
---|
1582 | if (!fp) { |
---|
1583 | fprintf(stderr,"Unable to open file input.130 in dirSample or Data/Sample directory\n"); |
---|
1584 | } else { |
---|
1585 | int problem; |
---|
1586 | char temp[100]; |
---|
1587 | // read and skip |
---|
1588 | fscanf(fp,"%s",temp); |
---|
1589 | assert (!strcmp(temp,"BEGIN")); |
---|
1590 | fscanf(fp,"%*s %*s %d %d %*s %*s %d %*s",&problem, &numberRows, |
---|
1591 | &numberColumns); |
---|
1592 | // scan down to SUPPLY |
---|
1593 | while (fgets(temp,100,fp)) { |
---|
1594 | if (!strncmp(temp,"SUPPLY",6)) |
---|
1595 | break; |
---|
1596 | } |
---|
1597 | if (strncmp(temp,"SUPPLY",6)) { |
---|
1598 | fprintf(stderr,"Unable to find SUPPLY\n"); |
---|
1599 | exit(2); |
---|
1600 | } |
---|
1601 | // get space for rhs |
---|
1602 | double * lower = new double[numberRows]; |
---|
1603 | double * upper = new double[numberRows]; |
---|
1604 | int i; |
---|
1605 | for (i=0;i<numberRows;i++) { |
---|
1606 | lower[i]=0.0; |
---|
1607 | upper[i]=0.0; |
---|
1608 | } |
---|
1609 | // ***** Remember to convert to C notation |
---|
1610 | while (fgets(temp,100,fp)) { |
---|
1611 | int row; |
---|
1612 | int value; |
---|
1613 | if (!strncmp(temp,"ARCS",4)) |
---|
1614 | break; |
---|
1615 | sscanf(temp,"%d %d",&row,&value); |
---|
1616 | upper[row-1]=-value; |
---|
1617 | lower[row-1]=-value; |
---|
1618 | } |
---|
1619 | if (strncmp(temp,"ARCS",4)) { |
---|
1620 | fprintf(stderr,"Unable to find ARCS\n"); |
---|
1621 | exit(2); |
---|
1622 | } |
---|
1623 | // number of columns may be underestimate |
---|
1624 | int * head = new int[2*numberColumns]; |
---|
1625 | int * tail = new int[2*numberColumns]; |
---|
1626 | int * ub = new int[2*numberColumns]; |
---|
1627 | int * cost = new int[2*numberColumns]; |
---|
1628 | // ***** Remember to convert to C notation |
---|
1629 | numberColumns=0; |
---|
1630 | while (fgets(temp,100,fp)) { |
---|
1631 | int iHead; |
---|
1632 | int iTail; |
---|
1633 | int iUb; |
---|
1634 | int iCost; |
---|
1635 | if (!strncmp(temp,"DEMAND",6)) |
---|
1636 | break; |
---|
1637 | sscanf(temp,"%d %d %d %d",&iHead,&iTail,&iCost,&iUb); |
---|
1638 | iHead--; |
---|
1639 | iTail--; |
---|
1640 | head[numberColumns]=iHead; |
---|
1641 | tail[numberColumns]=iTail; |
---|
1642 | ub[numberColumns]=iUb; |
---|
1643 | cost[numberColumns]=iCost; |
---|
1644 | numberColumns++; |
---|
1645 | } |
---|
1646 | if (strncmp(temp,"DEMAND",6)) { |
---|
1647 | fprintf(stderr,"Unable to find DEMAND\n"); |
---|
1648 | exit(2); |
---|
1649 | } |
---|
1650 | // ***** Remember to convert to C notation |
---|
1651 | while (fgets(temp,100,fp)) { |
---|
1652 | int row; |
---|
1653 | int value; |
---|
1654 | if (!strncmp(temp,"END",3)) |
---|
1655 | break; |
---|
1656 | sscanf(temp,"%d %d",&row,&value); |
---|
1657 | upper[row-1]=value; |
---|
1658 | lower[row-1]=value; |
---|
1659 | } |
---|
1660 | if (strncmp(temp,"END",3)) { |
---|
1661 | fprintf(stderr,"Unable to find END\n"); |
---|
1662 | exit(2); |
---|
1663 | } |
---|
1664 | printf("Problem %d has %d rows and %d columns\n",problem, |
---|
1665 | numberRows,numberColumns); |
---|
1666 | fclose(fp); |
---|
1667 | ClpSimplex model; |
---|
1668 | // now build model |
---|
1669 | |
---|
1670 | double * objective =new double[numberColumns]; |
---|
1671 | double * lowerColumn = new double[numberColumns]; |
---|
1672 | double * upperColumn = new double[numberColumns]; |
---|
1673 | |
---|
1674 | double * element = new double [2*numberColumns]; |
---|
1675 | CoinBigIndex * start = new CoinBigIndex [numberColumns+1]; |
---|
1676 | int * row = new int[2*numberColumns]; |
---|
1677 | start[numberColumns]=2*numberColumns; |
---|
1678 | for (i=0;i<numberColumns;i++) { |
---|
1679 | start[i]=2*i; |
---|
1680 | element[2*i]=-1.0; |
---|
1681 | element[2*i+1]=1.0; |
---|
1682 | row[2*i]=head[i]; |
---|
1683 | row[2*i+1]=tail[i]; |
---|
1684 | lowerColumn[i]=0.0; |
---|
1685 | upperColumn[i]=ub[i]; |
---|
1686 | objective[i]=cost[i]; |
---|
1687 | } |
---|
1688 | // Create Packed Matrix |
---|
1689 | CoinPackedMatrix matrix; |
---|
1690 | int * lengths = NULL; |
---|
1691 | matrix.assignMatrix(true,numberRows,numberColumns, |
---|
1692 | 2*numberColumns,element,row,start,lengths); |
---|
1693 | // load model |
---|
1694 | model.loadProblem(matrix, |
---|
1695 | lowerColumn,upperColumn,objective, |
---|
1696 | lower,upper); |
---|
1697 | model.factorization()->maximumPivots(200+model.numberRows()/100); |
---|
1698 | model.createStatus(); |
---|
1699 | double time1 = CoinCpuTime(); |
---|
1700 | model.dual(); |
---|
1701 | std::cout<<"Network problem, ClpPackedMatrix took "<<CoinCpuTime()-time1<<" seconds"<<std::endl; |
---|
1702 | ClpPlusMinusOneMatrix * plusMinus = new ClpPlusMinusOneMatrix(matrix); |
---|
1703 | assert (plusMinus->getIndices()); // would be zero if not +- one |
---|
1704 | //ClpPlusMinusOneMatrix *plusminus_matrix; |
---|
1705 | |
---|
1706 | //plusminus_matrix = new ClpPlusMinusOneMatrix; |
---|
1707 | |
---|
1708 | //plusminus_matrix->passInCopy(numberRows, numberColumns, true, plusMinus->getMutableIndices(), |
---|
1709 | // plusMinus->startPositive(),plusMinus->startNegative()); |
---|
1710 | model.loadProblem(*plusMinus, |
---|
1711 | lowerColumn,upperColumn,objective, |
---|
1712 | lower,upper); |
---|
1713 | //model.replaceMatrix( plusminus_matrix , true); |
---|
1714 | delete plusMinus; |
---|
1715 | //model.createStatus(); |
---|
1716 | //model.initialSolve(); |
---|
1717 | //model.writeMps("xx.mps"); |
---|
1718 | |
---|
1719 | model.factorization()->maximumPivots(200+model.numberRows()/100); |
---|
1720 | model.createStatus(); |
---|
1721 | time1 = CoinCpuTime(); |
---|
1722 | model.dual(); |
---|
1723 | std::cout<<"Network problem, ClpPlusMinusOneMatrix took "<<CoinCpuTime()-time1<<" seconds"<<std::endl; |
---|
1724 | ClpNetworkMatrix network(numberColumns,head,tail); |
---|
1725 | model.loadProblem(network, |
---|
1726 | lowerColumn,upperColumn,objective, |
---|
1727 | lower,upper); |
---|
1728 | |
---|
1729 | model.factorization()->maximumPivots(200+model.numberRows()/100); |
---|
1730 | model.createStatus(); |
---|
1731 | time1 = CoinCpuTime(); |
---|
1732 | model.dual(); |
---|
1733 | std::cout<<"Network problem, ClpNetworkMatrix took "<<CoinCpuTime()-time1<<" seconds"<<std::endl; |
---|
1734 | delete [] lower; |
---|
1735 | delete [] upper; |
---|
1736 | delete [] head; |
---|
1737 | delete [] tail; |
---|
1738 | delete [] ub; |
---|
1739 | delete [] cost; |
---|
1740 | delete [] objective; |
---|
1741 | delete [] lowerColumn; |
---|
1742 | delete [] upperColumn; |
---|
1743 | } |
---|
1744 | } |
---|
1745 | #ifdef QUADRATIC |
---|
1746 | // Test quadratic to solve linear |
---|
1747 | if (1) { |
---|
1748 | CoinMpsIO m; |
---|
1749 | std::string fn = dirSample+"exmip1"; |
---|
1750 | m.readMps(fn.c_str(),"mps"); |
---|
1751 | ClpSimplex solution; |
---|
1752 | solution.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(), |
---|
1753 | m.getObjCoefficients(), |
---|
1754 | m.getRowLower(),m.getRowUpper()); |
---|
1755 | //solution.dual(); |
---|
1756 | // get quadratic part |
---|
1757 | int numberColumns=solution.numberColumns(); |
---|
1758 | int * start=new int [numberColumns+1]; |
---|
1759 | int * column = new int[numberColumns]; |
---|
1760 | double * element = new double[numberColumns]; |
---|
1761 | int i; |
---|
1762 | start[0]=0; |
---|
1763 | int n=0; |
---|
1764 | int kk=numberColumns-1; |
---|
1765 | int kk2=numberColumns-1; |
---|
1766 | for (i=0;i<numberColumns;i++) { |
---|
1767 | if (i>=kk) { |
---|
1768 | column[n]=i; |
---|
1769 | if (i>=kk2) |
---|
1770 | element[n]=1.0e-1; |
---|
1771 | else |
---|
1772 | element[n]=0.0; |
---|
1773 | n++; |
---|
1774 | } |
---|
1775 | start[i+1]=n; |
---|
1776 | } |
---|
1777 | // Load up objective |
---|
1778 | solution.loadQuadraticObjective(numberColumns,start,column,element); |
---|
1779 | delete [] start; |
---|
1780 | delete [] column; |
---|
1781 | delete [] element; |
---|
1782 | //solution.quadraticSLP(50,1.0e-4); |
---|
1783 | CoinRelFltEq eq(1.0e-4); |
---|
1784 | //assert(eq(objValue,820.0)); |
---|
1785 | //solution.setLogLevel(63); |
---|
1786 | solution.primal(); |
---|
1787 | printSol(solution); |
---|
1788 | //assert(eq(objValue,3.2368421)); |
---|
1789 | //exit(77); |
---|
1790 | } |
---|
1791 | // Test quadratic |
---|
1792 | if (1) { |
---|
1793 | CoinMpsIO m; |
---|
1794 | std::string fn = dirSample+"share2qp"; |
---|
1795 | //fn = "share2qpb"; |
---|
1796 | m.readMps(fn.c_str(),"mps"); |
---|
1797 | ClpSimplex model; |
---|
1798 | model.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(), |
---|
1799 | m.getObjCoefficients(), |
---|
1800 | m.getRowLower(),m.getRowUpper()); |
---|
1801 | model.dual(); |
---|
1802 | // get quadratic part |
---|
1803 | int * start=NULL; |
---|
1804 | int * column = NULL; |
---|
1805 | double * element = NULL; |
---|
1806 | m.readQuadraticMps(NULL,start,column,element,2); |
---|
1807 | int column2[200]; |
---|
1808 | double element2[200]; |
---|
1809 | int start2[80]; |
---|
1810 | int j; |
---|
1811 | start2[0]=0; |
---|
1812 | int nel=0; |
---|
1813 | bool good=false; |
---|
1814 | for (j=0;j<79;j++) { |
---|
1815 | if (start[j]==start[j+1]) { |
---|
1816 | column2[nel]=j; |
---|
1817 | element2[nel]=0.0; |
---|
1818 | nel++; |
---|
1819 | } else { |
---|
1820 | int i; |
---|
1821 | for (i=start[j];i<start[j+1];i++) { |
---|
1822 | column2[nel]=column[i]; |
---|
1823 | element2[nel++]=element[i]; |
---|
1824 | } |
---|
1825 | } |
---|
1826 | start2[j+1]=nel; |
---|
1827 | } |
---|
1828 | // Load up objective |
---|
1829 | if (good) |
---|
1830 | model.loadQuadraticObjective(model.numberColumns(),start2,column2,element2); |
---|
1831 | else |
---|
1832 | model.loadQuadraticObjective(model.numberColumns(),start,column,element); |
---|
1833 | delete [] start; |
---|
1834 | delete [] column; |
---|
1835 | delete [] element; |
---|
1836 | int numberColumns=model.numberColumns(); |
---|
1837 | model.scaling(0); |
---|
1838 | #if 0 |
---|
1839 | model.nonlinearSLP(50,1.0e-4); |
---|
1840 | #else |
---|
1841 | // Get feasible |
---|
1842 | ClpObjective * saveObjective = model.objectiveAsObject()->clone(); |
---|
1843 | ClpLinearObjective zeroObjective(NULL,numberColumns); |
---|
1844 | model.setObjective(&zeroObjective); |
---|
1845 | model.dual(); |
---|
1846 | model.setObjective(saveObjective); |
---|
1847 | delete saveObjective; |
---|
1848 | #endif |
---|
1849 | //model.setLogLevel(63); |
---|
1850 | //exit(77); |
---|
1851 | model.setFactorizationFrequency(10); |
---|
1852 | model.primal(); |
---|
1853 | printSol(model); |
---|
1854 | #ifndef NDEBUG |
---|
1855 | double objValue = model.getObjValue(); |
---|
1856 | #endif |
---|
1857 | CoinRelFltEq eq(1.0e-4); |
---|
1858 | assert(eq(objValue,-400.92)); |
---|
1859 | // and again for barrier |
---|
1860 | model.barrier(false); |
---|
1861 | //printSol(model); |
---|
1862 | model.allSlackBasis(); |
---|
1863 | model.primal(); |
---|
1864 | //printSol(model); |
---|
1865 | } |
---|
1866 | if (0) { |
---|
1867 | CoinMpsIO m; |
---|
1868 | std::string fn = "./beale"; |
---|
1869 | //fn = "./jensen"; |
---|
1870 | m.readMps(fn.c_str(),"mps"); |
---|
1871 | ClpSimplex solution; |
---|
1872 | solution.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(), |
---|
1873 | m.getObjCoefficients(), |
---|
1874 | m.getRowLower(),m.getRowUpper()); |
---|
1875 | solution.setDblParam(ClpObjOffset,m.objectiveOffset()); |
---|
1876 | solution.dual(); |
---|
1877 | // get quadratic part |
---|
1878 | int * start=NULL; |
---|
1879 | int * column = NULL; |
---|
1880 | double * element = NULL; |
---|
1881 | m.readQuadraticMps(NULL,start,column,element,2); |
---|
1882 | // Load up objective |
---|
1883 | solution.loadQuadraticObjective(solution.numberColumns(),start,column,element); |
---|
1884 | delete [] start; |
---|
1885 | delete [] column; |
---|
1886 | delete [] element; |
---|
1887 | solution.primal(1); |
---|
1888 | solution.nonlinearSLP(50,1.0e-4); |
---|
1889 | double objValue = solution.getObjValue(); |
---|
1890 | CoinRelFltEq eq(1.0e-4); |
---|
1891 | assert(eq(objValue,0.5)); |
---|
1892 | solution.primal(); |
---|
1893 | objValue = solution.getObjValue(); |
---|
1894 | assert(eq(objValue,0.5)); |
---|
1895 | } |
---|
1896 | #endif |
---|
1897 | // Test CoinStructuredModel |
---|
1898 | { |
---|
1899 | |
---|
1900 | // Sub block |
---|
1901 | CoinModel sub; |
---|
1902 | { |
---|
1903 | // matrix data |
---|
1904 | //deliberate hiccup of 2 between 0 and 1 |
---|
1905 | CoinBigIndex start[5]={0,4,7,8,9}; |
---|
1906 | int length[5]={2,3,1,1,1}; |
---|
1907 | int rows[11]={0,2,-1,-1,0,1,2,0,1,2}; |
---|
1908 | double elements[11]={7.0,2.0,1.0e10,1.0e10,-2.0,1.0,-2.0,1,1,1}; |
---|
1909 | CoinPackedMatrix matrix(true,3,5,8,elements,rows,start,length); |
---|
1910 | // by row |
---|
1911 | matrix.reverseOrdering(); |
---|
1912 | const double * element = matrix.getElements(); |
---|
1913 | const int * column = matrix.getIndices(); |
---|
1914 | const CoinBigIndex * rowStart = matrix.getVectorStarts(); |
---|
1915 | const int * rowLength = matrix.getVectorLengths(); |
---|
1916 | |
---|
1917 | // rim data |
---|
1918 | //double objective[5]={-4.0,1.0,0.0,0.0,0.0}; |
---|
1919 | double rowLower[3]={14.0,3.0,3.0}; |
---|
1920 | double rowUpper[3]={14.0,3.0,3.0}; |
---|
1921 | //double columnLower[5]={0.0,0.0,0.0,0.0,0.0}; |
---|
1922 | //double columnUpper[5]={100.0,100.0,100.0,100.0,100.0}; |
---|
1923 | |
---|
1924 | for (int i=0;i<3;i++) { |
---|
1925 | sub.addRow(rowLength[i],column+rowStart[i], |
---|
1926 | element+rowStart[i],rowLower[i],rowUpper[i]); |
---|
1927 | } |
---|
1928 | //for (int i=0;i<5;i++) { |
---|
1929 | //sub.setColumnBounds(i,columnLower[i],columnUpper[i]); |
---|
1930 | //sub.setColumnObjective(i,objective[i]); |
---|
1931 | //} |
---|
1932 | sub.convertMatrix(); |
---|
1933 | } |
---|
1934 | // Top |
---|
1935 | CoinModel top; |
---|
1936 | { |
---|
1937 | // matrix data |
---|
1938 | CoinBigIndex start[5]={0,2,4,6,8}; |
---|
1939 | int length[5]={2,2,2,2,2}; |
---|
1940 | int rows[10]={0,1,0,1,0,1,0,1,0,1}; |
---|
1941 | double elements[10]={1,1,1,1,1,1,1,1,1,1}; |
---|
1942 | CoinPackedMatrix matrix(true,2,5,8,elements,rows,start,length); |
---|
1943 | // by row |
---|
1944 | matrix.reverseOrdering(); |
---|
1945 | const double * element = matrix.getElements(); |
---|
1946 | const int * column = matrix.getIndices(); |
---|
1947 | const CoinBigIndex * rowStart = matrix.getVectorStarts(); |
---|
1948 | const int * rowLength = matrix.getVectorLengths(); |
---|
1949 | |
---|
1950 | // rim data |
---|
1951 | double objective[5]={-4.0,1.0,0.0,0.0,0.0}; |
---|
1952 | //double rowLower[3]={14.0,3.0,3.0}; |
---|
1953 | //double rowUpper[3]={14.0,3.0,3.0}; |
---|
1954 | double columnLower[5]={0.0,0.0,0.0,0.0,0.0}; |
---|
1955 | double columnUpper[5]={100.0,100.0,100.0,100.0,100.0}; |
---|
1956 | |
---|
1957 | for (int i=0;i<2;i++) { |
---|
1958 | top.addRow(rowLength[i],column+rowStart[i], |
---|
1959 | element+rowStart[i], |
---|
1960 | -COIN_DBL_MAX,COIN_DBL_MAX); |
---|
1961 | } |
---|
1962 | for (int i=0;i<5;i++) { |
---|
1963 | top.setColumnBounds(i,columnLower[i],columnUpper[i]); |
---|
1964 | top.setColumnObjective(i,objective[i]); |
---|
1965 | } |
---|
1966 | top.convertMatrix(); |
---|
1967 | } |
---|
1968 | // Create a structured model |
---|
1969 | CoinStructuredModel structured; |
---|
1970 | int numberBlocks=5; |
---|
1971 | for (int i=0;i<numberBlocks;i++) { |
---|
1972 | std::string topName="row_master"; |
---|
1973 | std::string blockName="block_"; |
---|
1974 | char bName = 'a'+static_cast<char>(i); |
---|
1975 | blockName.append(1,bName); |
---|
1976 | structured.addBlock(topName,blockName,top); |
---|
1977 | structured.addBlock(blockName,blockName,sub); |
---|
1978 | } |
---|
1979 | // Set rhs on first block |
---|
1980 | CoinModel * first = structured.coinBlock(0); |
---|
1981 | for (int i=0;i<2;i++) { |
---|
1982 | first->setRowLower(i,0.0); |
---|
1983 | first->setRowUpper(i,100.0); |
---|
1984 | } |
---|
1985 | // Refresh whats set |
---|
1986 | structured.refresh(0); |
---|
1987 | // Could perturb stuff, but for first go don't bother |
---|
1988 | ClpSimplex fullModel; |
---|
1989 | // There is no original stuff set - think |
---|
1990 | fullModel.loadProblem(structured,false); |
---|
1991 | fullModel.dual(); |
---|
1992 | fullModel.dropNames(); |
---|
1993 | fullModel.writeMps("test.mps"); |
---|
1994 | // Make up very simple nested model - not realistic |
---|
1995 | // Create a structured model |
---|
1996 | CoinStructuredModel structured2; |
---|
1997 | numberBlocks=3; |
---|
1998 | for (int i=0;i<numberBlocks;i++) { |
---|
1999 | std::string blockName="block_"; |
---|
2000 | char bName = 'a'+static_cast<char>(i); |
---|
2001 | blockName.append(1,bName); |
---|
2002 | structured2.addBlock(blockName,blockName,structured); |
---|
2003 | } |
---|
2004 | fullModel.loadProblem(structured2,false); |
---|
2005 | fullModel.dual(); |
---|
2006 | fullModel.dropNames(); |
---|
2007 | fullModel.writeMps("test2.mps"); |
---|
2008 | } |
---|
2009 | } |
---|