source: trunk/Clp/src/unitTest.cpp @ 1525

Last change on this file since 1525 was 1525, checked in by mjs, 11 years ago

Formatted .cpp, .hpp, .c, .h files with "astyle -A4 -p". This matches the formatting used in the grand CBC reorganization.

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