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