source: branches/BSP/trunk/Cbc/src/unitTestClp.cpp @ 737

Last change on this file since 737 was 737, checked in by ladanyi, 14 years ago

Got rid of dependency on Data/... except for Sample

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