source: branches/devel/Cbc/src/unitTestClp.cpp @ 539

Last change on this file since 539 was 452, checked in by forrest, 13 years ago

modify CoinSolve? so will solve all 44 with -miplib
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 78.3 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        solveOptions=setupForSolve(iAlg,nameAlgorithm,0);
557        if (presolveOff)
558          solveOptions.setPresolveType(ClpSolve::presolveOff);
559      }
560      if (doVector) {
561        ClpMatrixBase * matrix = solution.clpMatrix();
562        if (dynamic_cast< ClpPackedMatrix*>(matrix)) {
563          ClpPackedMatrix * clpMatrix = dynamic_cast< ClpPackedMatrix*>(matrix);
564          clpMatrix->makeSpecialColumnCopy();
565        }
566      }
567      solution.initialSolve(solveOptions);
568      double time2 = CoinCpuTime()-time1;
569      timeTaken += time2;
570      printf("%s took %g seconds using algorithm %s\n",fn.c_str(),time2,nameAlgorithm.c_str());
571      // Test objective solution value
572      {
573        double soln = solution.objectiveValue();
574        CoinRelFltEq eq(objValueTol[m]);
575        std::cerr <<soln <<",  " <<objValue[m] <<" diff "<<
576          soln-objValue[m]<<std::endl;
577        if(!eq(soln,objValue[m]))
578          printf("** difference fails\n");
579      }
580    }
581    printf("Total time %g seconds\n",timeTaken);
582  }
583  else {
584    testingMessage( "***Skipped Testing on netlib    ***\n" );
585    testingMessage( "***use -netlib to test class***\n" );
586  }
587 
588  testingMessage( "All tests completed successfully\n" );
589  return 0;
590}
591
592 
593// Display message on stdout and stderr
594void testingMessage( const char * const msg )
595{
596  std::cerr <<msg;
597  //cout <<endl <<"*****************************************"
598  //     <<endl <<msg <<endl;
599}
600
601//--------------------------------------------------------------------------
602// test factorization methods and simplex method and simple barrier
603void
604ClpSimplexUnitTest(const std::string & mpsDir,
605                   const std::string & netlibDir)
606{
607 
608  CoinRelFltEq eq(0.000001);
609
610  {
611    ClpSimplex solution;
612 
613    // matrix data
614    //deliberate hiccup of 2 between 0 and 1
615    CoinBigIndex start[5]={0,4,7,8,9};
616    int length[5]={2,3,1,1,1};
617    int rows[11]={0,2,-1,-1,0,1,2,0,1,2};
618    double elements[11]={7.0,2.0,1.0e10,1.0e10,-2.0,1.0,-2.0,1,1,1};
619    CoinPackedMatrix matrix(true,3,5,8,elements,rows,start,length);
620   
621    // rim data
622    double objective[7]={-4.0,1.0,0.0,0.0,0.0,0.0,0.0};
623    double rowLower[5]={14.0,3.0,3.0,1.0e10,1.0e10};
624    double rowUpper[5]={14.0,3.0,3.0,-1.0e10,-1.0e10};
625    double colLower[7]={0.0,0.0,0.0,0.0,0.0,0.0,0.0};
626    double colUpper[7]={100.0,100.0,100.0,100.0,100.0,100.0,100.0};
627   
628    // basis 1
629    int rowBasis1[3]={-1,-1,-1};
630    int colBasis1[5]={1,1,-1,-1,1};
631    solution.loadProblem(matrix,colLower,colUpper,objective,
632                         rowLower,rowUpper);
633    int i;
634    solution.createStatus();
635    for (i=0;i<3;i++) {
636      if (rowBasis1[i]<0) {
637        solution.setRowStatus(i,ClpSimplex::atLowerBound);
638      } else {
639        solution.setRowStatus(i,ClpSimplex::basic);
640      }
641    }
642    for (i=0;i<5;i++) {
643      if (colBasis1[i]<0) {
644        solution.setColumnStatus(i,ClpSimplex::atLowerBound);
645      } else {
646        solution.setColumnStatus(i,ClpSimplex::basic);
647      }
648    }
649    solution.setLogLevel(3+4+8+16+32);
650    solution.primal();
651    for (i=0;i<3;i++) {
652      if (rowBasis1[i]<0) {
653        solution.setRowStatus(i,ClpSimplex::atLowerBound);
654      } else {
655        solution.setRowStatus(i,ClpSimplex::basic);
656      }
657    }
658    for (i=0;i<5;i++) {
659      if (colBasis1[i]<0) {
660        solution.setColumnStatus(i,ClpSimplex::atLowerBound);
661      } else {
662        solution.setColumnStatus(i,ClpSimplex::basic);
663      }
664    }
665    // intricate stuff does not work with scaling
666    solution.scaling(0);
667    int returnCode = solution.factorize ( );
668    assert(!returnCode);
669    const double * colsol = solution.primalColumnSolution();
670    const double * rowsol = solution.primalRowSolution();
671    solution.getSolution(rowsol,colsol);
672    double colsol1[5]={20.0/7.0,3.0,0.0,0.0,23.0/7.0};
673    for (i=0;i<5;i++) {
674      assert(eq(colsol[i],colsol1[i]));
675    }
676    // now feed in again without actually doing factorization
677    ClpFactorization factorization2 = *solution.factorization();
678    ClpSimplex solution2 = solution;
679    solution2.setFactorization(factorization2);
680    solution2.createStatus();
681    for (i=0;i<3;i++) {
682      if (rowBasis1[i]<0) {
683        solution2.setRowStatus(i,ClpSimplex::atLowerBound);
684      } else {
685        solution2.setRowStatus(i,ClpSimplex::basic);
686      }
687    }
688    for (i=0;i<5;i++) {
689      if (colBasis1[i]<0) {
690        solution2.setColumnStatus(i,ClpSimplex::atLowerBound);
691      } else {
692        solution2.setColumnStatus(i,ClpSimplex::basic);
693      }
694    }
695    // intricate stuff does not work with scaling
696    solution2.scaling(0);
697    solution2.getSolution(rowsol,colsol);
698    colsol = solution2.primalColumnSolution();
699    rowsol = solution2.primalRowSolution();
700    for (i=0;i<5;i++) {
701      assert(eq(colsol[i],colsol1[i]));
702    }
703    solution2.setDualBound(0.1);
704    solution2.dual();
705    objective[2]=-1.0;
706    objective[3]=-0.5;
707    objective[4]=10.0;
708    solution.dual();
709    for (i=0;i<3;i++) {
710      rowLower[i]=-1.0e20;
711      colUpper[i+2]=0.0;
712    }
713    solution.setLogLevel(3);
714    solution.dual();
715    double rowObjective[]={1.0,0.5,-10.0};
716    solution.loadProblem(matrix,colLower,colUpper,objective,
717                         rowLower,rowUpper,rowObjective);
718    solution.dual();
719    solution.loadProblem(matrix,colLower,colUpper,objective,
720                         rowLower,rowUpper,rowObjective);
721    solution.primal();
722  }
723#ifndef COIN_NO_CLP_MESSAGE
724  {   
725    CoinMpsIO m;
726    std::string fn = mpsDir+"exmip1";
727    m.readMps(fn.c_str(),"mps");
728    ClpSimplex solution;
729    solution.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
730                         m.getObjCoefficients(),
731                         m.getRowLower(),m.getRowUpper());
732    solution.dual();
733    // Test event handling
734    MyEventHandler handler;
735    solution.passInEventHandler(&handler);
736    int numberRows=solution.numberRows();
737    // make sure values pass has something to do
738    for (int i=0;i<numberRows;i++)
739      solution.setRowStatus(i,ClpSimplex::basic);
740    solution.primal(1);
741    assert (solution.secondaryStatus()==102); // Came out at end of pass
742  }
743  // Test Message handler
744  {   
745    CoinMpsIO m;
746    std::string fn = mpsDir+"exmip1";
747    //fn = "Test/subGams4";
748    m.readMps(fn.c_str(),"mps");
749    ClpSimplex model;
750    model.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
751                         m.getObjCoefficients(),
752                         m.getRowLower(),m.getRowUpper());
753    // Message handler
754    MyMessageHandler messageHandler(&model);
755    std::cout<<"Testing derived message handler"<<std::endl;
756    model.passInMessageHandler(&messageHandler);
757    model.messagesPointer()->setDetailMessage(1,102);
758    model.setFactorizationFrequency(10);
759    model.primal();
760
761    // Write saved solutions
762    int nc = model.getNumCols();
763    int s; 
764    std::deque<StdVectorDouble> fep = messageHandler.getFeasibleExtremePoints();
765    int numSavedSolutions = fep.size();
766    for ( s=0; s<numSavedSolutions; ++s ) {
767      const StdVectorDouble & solnVec = fep[s];
768      for ( int c=0; c<nc; ++c ) {
769        if (fabs(solnVec[c])>1.0e-8)
770          std::cout <<"Saved Solution: " <<s <<" ColNum: " <<c <<" Value: " <<solnVec[c] <<std::endl;
771      }
772    }
773    // Solve again without scaling
774    // and maximize then minimize
775    messageHandler.clearFeasibleExtremePoints();
776    model.scaling(0);
777    model.setOptimizationDirection(-1);
778    model.primal();
779    model.setOptimizationDirection(1);
780    model.primal();
781    fep = messageHandler.getFeasibleExtremePoints();
782    numSavedSolutions = fep.size();
783    for ( s=0; s<numSavedSolutions; ++s ) {
784      const StdVectorDouble & solnVec = fep[s];
785      for ( int c=0; c<nc; ++c ) {
786        if (fabs(solnVec[c])>1.0e-8)
787          std::cout <<"Saved Solution: " <<s <<" ColNum: " <<c <<" Value: " <<solnVec[c] <<std::endl;
788      }
789    }
790  }
791#endif
792  // Test dual ranging
793  {   
794    CoinMpsIO m;
795    std::string fn = mpsDir+"exmip1";
796    m.readMps(fn.c_str(),"mps");
797    ClpSimplex model;
798    model.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
799                         m.getObjCoefficients(),
800                         m.getRowLower(),m.getRowUpper());
801    model.primal();
802    int which[13] = {0,1,2,3,4,5,6,7,8,9,10,11,12};
803    double costIncrease[13];
804    int sequenceIncrease[13];
805    double costDecrease[13];
806    int sequenceDecrease[13];
807    // ranging
808    model.dualRanging(13,which,costIncrease,sequenceIncrease,
809                      costDecrease,sequenceDecrease);
810    int i;
811    for ( i=0;i<13;i++)
812      printf("%d increase %g %d, decrease %g %d\n",
813             i,costIncrease[i],sequenceIncrease[i],
814             costDecrease[i],sequenceDecrease[i]);
815    assert (fabs(costDecrease[3])<1.0e-4);
816    assert (fabs(costIncrease[7]-1.0)<1.0e-4);
817    model.setOptimizationDirection(-1);
818    {
819      int j;
820      double * obj = model.objective();
821      int n=model.numberColumns();
822      for (j=0;j<n;j++) 
823        obj[j] *= -1.0;
824    }
825    double costIncrease2[13];
826    int sequenceIncrease2[13];
827    double costDecrease2[13];
828    int sequenceDecrease2[13];
829    // ranging
830    model.dualRanging(13,which,costIncrease2,sequenceIncrease2,
831                      costDecrease2,sequenceDecrease2);
832    for (i=0;i<13;i++) {
833      assert (fabs(costIncrease[i]-costDecrease2[i])<1.0e-6);
834      assert (fabs(costDecrease[i]-costIncrease2[i])<1.0e-6);
835      assert (sequenceIncrease[i]==sequenceDecrease2[i]);
836      assert (sequenceDecrease[i]==sequenceIncrease2[i]);
837    }
838    // Now delete all rows and see what happens
839    model.deleteRows(model.numberRows(),which);
840    model.primal();
841    // ranging
842    if (!model.dualRanging(8,which,costIncrease,sequenceIncrease,
843                           costDecrease,sequenceDecrease)) {
844      for (i=0;i<8;i++) {
845        printf("%d increase %g %d, decrease %g %d\n",
846               i,costIncrease[i],sequenceIncrease[i],
847               costDecrease[i],sequenceDecrease[i]);
848      }
849    }
850  }
851  // Test primal ranging
852  {   
853    CoinMpsIO m;
854    std::string fn = mpsDir+"exmip1";
855    m.readMps(fn.c_str(),"mps");
856    ClpSimplex model;
857    model.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
858                         m.getObjCoefficients(),
859                         m.getRowLower(),m.getRowUpper());
860    model.primal();
861    int which[13] = {0,1,2,3,4,5,6,7,8,9,10,11,12};
862    double valueIncrease[13];
863    int sequenceIncrease[13];
864    double valueDecrease[13];
865    int sequenceDecrease[13];
866    // ranging
867    model.primalRanging(13,which,valueIncrease,sequenceIncrease,
868                      valueDecrease,sequenceDecrease);
869    int i;
870    for ( i=0;i<13;i++)
871      printf("%d increase %g %d, decrease %g %d\n",
872             i,valueIncrease[i],sequenceIncrease[i],
873             valueDecrease[i],sequenceDecrease[i]);
874    assert (fabs(valueDecrease[3]-0.642857)<1.0e-4);
875    assert (fabs(valueDecrease[8]-2.95113)<1.0e-4);
876#if 0
877    // out until I find optimization bug
878    // Test parametrics
879    ClpSimplexOther * model2 = (ClpSimplexOther *) (&model);
880    double rhs[]={ 1.0,2.0,3.0,4.0,5.0};
881    double endingTheta=1.0;
882    model2->scaling(0);
883    model2->setLogLevel(63);
884    model2->parametrics(0.0,endingTheta,0.1,
885                        NULL,NULL,rhs,rhs,NULL);
886#endif
887  }
888  // Test binv etc
889  {   
890    /*
891       Wolsey : Page 130
892       max 4x1 -  x2
893       7x1 - 2x2    <= 14
894       x2    <= 3
895       2x1 - 2x2    <= 3
896       x1 in Z+, x2 >= 0
897
898       note slacks are -1 in Clp so signs may be different
899    */
900   
901    int n_cols = 2;
902    int n_rows = 3;
903   
904    double obj[2] = {-4.0, 1.0};
905    double collb[2] = {0.0, 0.0};
906    double colub[2] = {COIN_DBL_MAX, COIN_DBL_MAX};
907    double rowlb[3] = {-COIN_DBL_MAX, -COIN_DBL_MAX, -COIN_DBL_MAX};
908    double rowub[3] = {14.0, 3.0, 3.0};
909   
910    int rowIndices[5] =  {0,     2,    0,    1,    2};
911    int colIndices[5] =  {0,     0,    1,    1,    1};
912    double elements[5] = {7.0, 2.0, -2.0,  1.0, -2.0};
913    CoinPackedMatrix M(true, rowIndices, colIndices, elements, 5);
914
915    ClpSimplex model;
916    model.loadProblem(M, collb, colub, obj, rowlb, rowub);
917    model.dual(0,1); // keep factorization
918   
919    //check that the tableau matches wolsey (B-1 A)
920    // slacks in second part of binvA
921    double * binvA = (double*) malloc((n_cols+n_rows) * sizeof(double));
922   
923    printf("B-1 A by row\n");
924    int i;
925    for( i = 0; i < n_rows; i++){
926      model.getBInvARow(i, binvA,binvA+n_cols);
927      printf("row: %d -> ",i);
928      for(int j=0; j < n_cols+n_rows; j++){
929        printf("%g, ", binvA[j]);
930      }
931      printf("\n");
932    }
933    // See if can re-use factorization AND arrays
934    model.primal(0,3+4); // keep factorization
935    // And do by column
936    printf("B-1 A by column\n");
937    for( i = 0; i < n_rows+n_cols; i++){
938      model.getBInvACol(i, binvA);
939      printf("column: %d -> ",i);
940      for(int j=0; j < n_rows; j++){
941        printf("%g, ", binvA[j]);
942      }
943      printf("\n");
944    }
945    free(binvA);
946    model.setColUpper(1,2.0);
947    model.dual(0,2+4); // use factorization and arrays
948    model.dual(0,2); // hopefully will not use factorization
949    model.primal(0,3+4); // keep factorization
950    // but say basis has changed
951    model.setWhatsChanged(model.whatsChanged()&(~512));
952    model.dual(0,2); // hopefully will not use factorization
953  }
954  // test steepest edge
955  {   
956    CoinMpsIO m;
957    std::string fn = netlibDir+"finnis";
958    int returnCode = m.readMps(fn.c_str(),"mps");
959    if (returnCode) {
960      // probable cause is that gz not there
961      fprintf(stderr,"Unable to open finnis.mps in Data/Netlib!\n");
962      fprintf(stderr,"Most probable cause is finnis.mps is gzipped i.e. finnis.mps.gz and libz has not been activated\n");
963      fprintf(stderr,"Either gunzip files or edit Makefiles/Makefile.location to get libz\n");
964      exit(999);
965    }
966    ClpModel model;
967    model.loadProblem(*m.getMatrixByCol(),m.getColLower(),
968                    m.getColUpper(),
969                    m.getObjCoefficients(),
970                    m.getRowLower(),m.getRowUpper());
971    ClpSimplex solution(model);
972
973    solution.scaling(1); 
974    solution.setDualBound(1.0e8);
975    //solution.factorization()->maximumPivots(1);
976    //solution.setLogLevel(3);
977    solution.setDualTolerance(1.0e-7);
978    // set objective sense,
979    ClpDualRowSteepest steep;
980    solution.setDualRowPivotAlgorithm(steep);
981    solution.setDblParam(ClpObjOffset,m.objectiveOffset());
982    solution.dual();
983  }
984  // test normal solution
985  {   
986    CoinMpsIO m;
987    std::string fn = netlibDir+"afiro";
988    m.readMps(fn.c_str(),"mps");
989    ClpSimplex solution;
990    ClpModel model;
991    // do twice - without and with scaling
992    int iPass;
993    for (iPass=0;iPass<2;iPass++) {
994      // explicit row objective for testing
995      int nr = m.getNumRows();
996      double * rowObj = new double[nr];
997      CoinFillN(rowObj,nr,0.0);
998      model.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
999                      m.getObjCoefficients(),
1000                      m.getRowLower(),m.getRowUpper(),rowObj);
1001      delete [] rowObj;
1002      solution = ClpSimplex(model);
1003      if (iPass) {
1004        solution.scaling();
1005      }
1006      solution.dual();
1007      solution.dual();
1008      // test optimal
1009      assert (solution.status()==0);
1010      int numberColumns = solution.numberColumns();
1011      int numberRows = solution.numberRows();
1012      CoinPackedVector colsol(numberColumns,solution.primalColumnSolution());
1013      double * objective = solution.objective();
1014      double objValue = colsol.dotProduct(objective);
1015      CoinRelFltEq eq(1.0e-8);
1016      assert(eq(objValue,-4.6475314286e+02));
1017      // Test auxiliary model
1018      //solution.scaling(0);
1019      solution.auxiliaryModel(63-2); // bounds may change
1020      solution.dual();
1021      solution.primal();
1022      solution.allSlackBasis();
1023      solution.dual();
1024      assert(eq(solution.objectiveValue(),-4.6475314286e+02));
1025      solution.auxiliaryModel(-1);
1026      solution.dual();
1027      assert(eq(solution.objectiveValue(),-4.6475314286e+02));
1028      double * lower = solution.columnLower();
1029      double * upper = solution.columnUpper();
1030      double * sol = solution.primalColumnSolution();
1031      double * result = new double[numberColumns];
1032      CoinFillN ( result, numberColumns,0.0);
1033      solution.matrix()->transposeTimes(solution.dualRowSolution(), result);
1034      int iRow , iColumn;
1035      // see if feasible and dual feasible
1036      for (iColumn=0;iColumn<numberColumns;iColumn++) {
1037        double value = sol[iColumn];
1038        assert(value<upper[iColumn]+1.0e-8);
1039        assert(value>lower[iColumn]-1.0e-8);
1040        value = objective[iColumn]-result[iColumn];
1041        assert (value>-1.0e-5);
1042        if (sol[iColumn]>1.0e-5)
1043          assert (value<1.0e-5);
1044      }
1045      delete [] result;
1046      result = new double[numberRows];
1047      CoinFillN ( result, numberRows,0.0);
1048      solution.matrix()->times(colsol, result);
1049      lower = solution.rowLower();
1050      upper = solution.rowUpper();
1051      sol = solution.primalRowSolution();
1052      for (iRow=0;iRow<numberRows;iRow++) {
1053        double value = result[iRow];
1054        assert(eq(value,sol[iRow]));
1055        assert(value<upper[iRow]+1.0e-8);
1056        assert(value>lower[iRow]-1.0e-8);
1057      }
1058      delete [] result;
1059      // test row objective
1060      double * rowObjective = solution.rowObjective();
1061      CoinDisjointCopyN(solution.dualRowSolution(),numberRows,rowObjective);
1062      CoinDisjointCopyN(solution.dualColumnSolution(),numberColumns,objective);
1063      // this sets up all slack basis
1064      solution.createStatus();
1065      solution.dual();
1066      CoinFillN(rowObjective,numberRows,0.0);
1067      CoinDisjointCopyN(m.getObjCoefficients(),numberColumns,objective);
1068      solution.dual();
1069    }
1070  }
1071  // test unbounded
1072  {   
1073    CoinMpsIO m;
1074    std::string fn = netlibDir+"brandy";
1075    m.readMps(fn.c_str(),"mps");
1076    ClpSimplex solution;
1077    // do twice - without and with scaling
1078    int iPass;
1079    for (iPass=0;iPass<2;iPass++) {
1080      solution.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
1081                      m.getObjCoefficients(),
1082                      m.getRowLower(),m.getRowUpper());
1083      if (iPass)
1084        solution.scaling();
1085      solution.setOptimizationDirection(-1);
1086      // test unbounded and ray
1087#ifdef DUAL
1088      solution.setDualBound(100.0);
1089      solution.dual();
1090#else
1091      solution.primal();
1092#endif
1093      assert (solution.status()==2);
1094      int numberColumns = solution.numberColumns();
1095      int numberRows = solution.numberRows();
1096      double * lower = solution.columnLower();
1097      double * upper = solution.columnUpper();
1098      double * sol = solution.primalColumnSolution();
1099      double * ray = solution.unboundedRay();
1100      double * objective = solution.objective();
1101      double objChange=0.0;
1102      int iRow , iColumn;
1103      // make sure feasible and columns form ray
1104      for (iColumn=0;iColumn<numberColumns;iColumn++) {
1105        double value = sol[iColumn];
1106        assert(value<upper[iColumn]+1.0e-8);
1107        assert(value>lower[iColumn]-1.0e-8);
1108        value = ray[iColumn];
1109        if (value>0.0)
1110          assert(upper[iColumn]>1.0e30);
1111        else if (value<0.0)
1112          assert(lower[iColumn]<-1.0e30);
1113        objChange += value*objective[iColumn];
1114      }
1115      // make sure increasing objective
1116      assert(objChange>0.0);
1117      double * result = new double[numberRows];
1118      CoinFillN ( result, numberRows,0.0);
1119      solution.matrix()->times(sol, result);
1120      lower = solution.rowLower();
1121      upper = solution.rowUpper();
1122      sol = solution.primalRowSolution();
1123      for (iRow=0;iRow<numberRows;iRow++) {
1124        double value = result[iRow];
1125        assert(eq(value,sol[iRow]));
1126        assert(value<upper[iRow]+2.0e-8);
1127        assert(value>lower[iRow]-2.0e-8);
1128      }
1129      CoinFillN ( result, numberRows,0.0);
1130      solution.matrix()->times(ray, result);
1131      // there may be small differences (especially if scaled)
1132      for (iRow=0;iRow<numberRows;iRow++) {
1133        double value = result[iRow];
1134        if (value>1.0e-8)
1135          assert(upper[iRow]>1.0e30);
1136        else if (value<-1.0e-8)
1137          assert(lower[iRow]<-1.0e30);
1138      }
1139      delete [] result;
1140      delete [] ray;
1141    }
1142  }
1143  // test infeasible
1144  {   
1145    CoinMpsIO m;
1146    std::string fn = netlibDir+"brandy";
1147    m.readMps(fn.c_str(),"mps");
1148    ClpSimplex solution;
1149    // do twice - without and with scaling
1150    int iPass;
1151    for (iPass=0;iPass<2;iPass++) {
1152      solution.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
1153                      m.getObjCoefficients(),
1154                      m.getRowLower(),m.getRowUpper());
1155      if (iPass)
1156        solution.scaling();
1157      // test infeasible and ray
1158      solution.columnUpper()[0]=0.0;
1159#ifdef DUAL
1160      solution.setDualBound(100.0);
1161      solution.dual();
1162#else
1163      solution.primal();
1164#endif
1165      assert (solution.status()==1);
1166      int numberColumns = solution.numberColumns();
1167      int numberRows = solution.numberRows();
1168      double * lower = solution.rowLower();
1169      double * upper = solution.rowUpper();
1170      double * ray = solution.infeasibilityRay();
1171      assert(ray);
1172      // construct proof of infeasibility
1173      int iRow , iColumn;
1174      double lo=0.0,up=0.0;
1175      int nl=0,nu=0;
1176      for (iRow=0;iRow<numberRows;iRow++) {
1177        if (lower[iRow]>-1.0e20) {
1178          lo += ray[iRow]*lower[iRow];
1179        } else {
1180          if (ray[iRow]>1.0e-8) 
1181            nl++;
1182        }
1183        if (upper[iRow]<1.0e20) {
1184          up += ray[iRow]*upper[iRow];
1185        } else {
1186          if (ray[iRow]>1.0e-8) 
1187            nu++;
1188        }
1189      }
1190      if (nl)
1191        lo=-1.0e100;
1192      if (nu)
1193        up=1.0e100;
1194      double * result = new double[numberColumns];
1195      double lo2=0.0,up2=0.0;
1196      CoinFillN ( result, numberColumns,0.0);
1197      solution.matrix()->transposeTimes(ray, result);
1198      lower = solution.columnLower();
1199      upper = solution.columnUpper();
1200      nl=nu=0;
1201      for (iColumn=0;iColumn<numberColumns;iColumn++) {
1202        if (result[iColumn]>1.0e-8) {
1203          if (lower[iColumn]>-1.0e20)
1204            lo2 += result[iColumn]*lower[iColumn];
1205          else
1206            nl++;
1207          if (upper[iColumn]<1.0e20)
1208            up2 += result[iColumn]*upper[iColumn];
1209          else
1210            nu++;
1211        } else if (result[iColumn]<-1.0e-8) {
1212          if (lower[iColumn]>-1.0e20)
1213            up2 += result[iColumn]*lower[iColumn];
1214          else
1215            nu++;
1216          if (upper[iColumn]<1.0e20)
1217            lo2 += result[iColumn]*upper[iColumn];
1218          else
1219            nl++;
1220        }
1221      }
1222      if (nl)
1223        lo2=-1.0e100;
1224      if (nu)
1225        up2=1.0e100;
1226      // make sure inconsistency
1227      assert(lo2>up||up2<lo);
1228      delete [] result;
1229      delete [] ray;
1230    }
1231  }
1232  // test delete and add
1233  {   
1234    CoinMpsIO m;
1235    std::string fn = netlibDir+"brandy";
1236    m.readMps(fn.c_str(),"mps");
1237    ClpSimplex solution;
1238    solution.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
1239                         m.getObjCoefficients(),
1240                         m.getRowLower(),m.getRowUpper());
1241    solution.dual();
1242    CoinRelFltEq eq(1.0e-8);
1243    assert(eq(solution.objectiveValue(),1.5185098965e+03));
1244
1245    int numberColumns = solution.numberColumns();
1246    int numberRows = solution.numberRows();
1247    double * saveObj = new double [numberColumns];
1248    double * saveLower = new double[numberRows+numberColumns];
1249    double * saveUpper = new double[numberRows+numberColumns];
1250    int * which = new int [numberRows+numberColumns];
1251
1252    int numberElements = m.getMatrixByCol()->getNumElements();
1253    int * starts = new int[numberRows+numberColumns];
1254    int * index = new int[numberElements];
1255    double * element = new double[numberElements];
1256
1257    const CoinBigIndex * startM;
1258    const int * lengthM;
1259    const int * indexM;
1260    const double * elementM;
1261
1262    int n,nel;
1263
1264    // delete non basic columns
1265    n=0;
1266    nel=0;
1267    int iRow , iColumn;
1268    const double * lower = m.getColLower();
1269    const double * upper = m.getColUpper();
1270    const double * objective = m.getObjCoefficients();
1271    startM = m.getMatrixByCol()->getVectorStarts();
1272    lengthM = m.getMatrixByCol()->getVectorLengths();
1273    indexM = m.getMatrixByCol()->getIndices();
1274    elementM = m.getMatrixByCol()->getElements();
1275    starts[0]=0;
1276    for (iColumn=0;iColumn<numberColumns;iColumn++) {
1277      if (solution.getColumnStatus(iColumn)!=ClpSimplex::basic) {
1278        saveObj[n]=objective[iColumn];
1279        saveLower[n]=lower[iColumn];
1280        saveUpper[n]=upper[iColumn];
1281        int j;
1282        for (j=startM[iColumn];j<startM[iColumn]+lengthM[iColumn];j++) {
1283          index[nel]=indexM[j];
1284          element[nel++]=elementM[j];
1285        }
1286        which[n++]=iColumn;
1287        starts[n]=nel;
1288      }
1289    }
1290    solution.deleteColumns(n,which);
1291    solution.dual();
1292    // Put back
1293    solution.addColumns(n,saveLower,saveUpper,saveObj,
1294                        starts,index,element);
1295    solution.dual();
1296    assert(eq(solution.objectiveValue(),1.5185098965e+03));
1297    // Delete all columns and add back
1298    n=0;
1299    nel=0;
1300    starts[0]=0;
1301    lower = m.getColLower();
1302    upper = m.getColUpper();
1303    objective = m.getObjCoefficients();
1304    for (iColumn=0;iColumn<numberColumns;iColumn++) {
1305      saveObj[n]=objective[iColumn];
1306      saveLower[n]=lower[iColumn];
1307      saveUpper[n]=upper[iColumn];
1308      int j;
1309      for (j=startM[iColumn];j<startM[iColumn]+lengthM[iColumn];j++) {
1310        index[nel]=indexM[j];
1311        element[nel++]=elementM[j];
1312      }
1313      which[n++]=iColumn;
1314      starts[n]=nel;
1315    }
1316    solution.deleteColumns(n,which);
1317    solution.dual();
1318    // Put back
1319    solution.addColumns(n,saveLower,saveUpper,saveObj,
1320                        starts,index,element);
1321    solution.dual();
1322    assert(eq(solution.objectiveValue(),1.5185098965e+03));
1323
1324    // reload with original
1325    solution.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
1326                         m.getObjCoefficients(),
1327                         m.getRowLower(),m.getRowUpper());
1328    // delete half rows
1329    n=0;
1330    nel=0;
1331    lower = m.getRowLower();
1332    upper = m.getRowUpper();
1333    startM = m.getMatrixByRow()->getVectorStarts();
1334    lengthM = m.getMatrixByRow()->getVectorLengths();
1335    indexM = m.getMatrixByRow()->getIndices();
1336    elementM = m.getMatrixByRow()->getElements();
1337    starts[0]=0;
1338    for (iRow=0;iRow<numberRows;iRow++) {
1339      if ((iRow&1)==0) {
1340        saveLower[n]=lower[iRow];
1341        saveUpper[n]=upper[iRow];
1342        int j;
1343        for (j=startM[iRow];j<startM[iRow]+lengthM[iRow];j++) {
1344          index[nel]=indexM[j];
1345          element[nel++]=elementM[j];
1346        }
1347        which[n++]=iRow;
1348        starts[n]=nel;
1349      }
1350    }
1351    solution.deleteRows(n,which);
1352    solution.dual();
1353    // Put back
1354    solution.addRows(n,saveLower,saveUpper,
1355                        starts,index,element);
1356    solution.dual();
1357    assert(eq(solution.objectiveValue(),1.5185098965e+03));
1358    solution.writeMps("yy.mps");
1359    // Delete all rows
1360    n=0;
1361    nel=0;
1362    lower = m.getRowLower();
1363    upper = m.getRowUpper();
1364    starts[0]=0;
1365    for (iRow=0;iRow<numberRows;iRow++) {
1366      saveLower[n]=lower[iRow];
1367      saveUpper[n]=upper[iRow];
1368      int j;
1369      for (j=startM[iRow];j<startM[iRow]+lengthM[iRow];j++) {
1370        index[nel]=indexM[j];
1371        element[nel++]=elementM[j];
1372      }
1373      which[n++]=iRow;
1374      starts[n]=nel;
1375    }
1376    solution.deleteRows(n,which);
1377    solution.dual();
1378    // Put back
1379    solution.addRows(n,saveLower,saveUpper,
1380                        starts,index,element);
1381    solution.dual();
1382    solution.writeMps("xx.mps");
1383    assert(eq(solution.objectiveValue(),1.5185098965e+03));
1384    // Zero out status array to give some interest
1385    memset(solution.statusArray()+numberColumns,0,numberRows);
1386    solution.primal(1);
1387    assert(eq(solution.objectiveValue(),1.5185098965e+03));
1388    // Delete all columns and rows
1389    n=0;
1390    for (iColumn=0;iColumn<numberColumns;iColumn++) {
1391      which[n++]=iColumn;
1392      starts[n]=nel;
1393    }
1394    solution.deleteColumns(n,which);
1395    n=0;
1396    for (iRow=0;iRow<numberRows;iRow++) {
1397      which[n++]=iRow;
1398      starts[n]=nel;
1399    }
1400    solution.deleteRows(n,which);
1401
1402    delete [] saveObj;
1403    delete [] saveLower;
1404    delete [] saveUpper;
1405    delete [] which;
1406    delete [] starts;
1407    delete [] index;
1408    delete [] element;
1409  }
1410#if 1
1411  // Test barrier
1412  {
1413    CoinMpsIO m;
1414    std::string fn = mpsDir+"exmip1";
1415    m.readMps(fn.c_str(),"mps");
1416    ClpInterior solution;
1417    solution.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
1418                         m.getObjCoefficients(),
1419                         m.getRowLower(),m.getRowUpper());
1420    solution.primalDual();
1421  }
1422#endif
1423  // test network
1424#define QUADRATIC
1425  if (1) {   
1426    std::string fn = mpsDir+"input.130";
1427    int numberColumns;
1428    int numberRows;
1429   
1430    FILE * fp = fopen(fn.c_str(),"r");
1431    if (!fp) {
1432      // Try in Data
1433      fn = "Data/Sample/input.130";
1434      fp = fopen(fn.c_str(),"r");
1435    }
1436    if (!fp) {
1437      fprintf(stderr,"Unable to open file input.130 in mpsDir or Data/Sample directory\n");
1438    } else {
1439      int problem;
1440      char temp[100];
1441      // read and skip
1442      fscanf(fp,"%s",temp);
1443      assert (!strcmp(temp,"BEGIN"));
1444      fscanf(fp,"%*s %*s %d %d %*s %*s %d %*s",&problem, &numberRows, 
1445             &numberColumns);
1446      // scan down to SUPPLY
1447      while (fgets(temp,100,fp)) {
1448        if (!strncmp(temp,"SUPPLY",6))
1449          break;
1450      }
1451      if (strncmp(temp,"SUPPLY",6)) {
1452        fprintf(stderr,"Unable to find SUPPLY\n");
1453        exit(2);
1454      }
1455      // get space for rhs
1456      double * lower = new double[numberRows];
1457      double * upper = new double[numberRows];
1458      int i;
1459      for (i=0;i<numberRows;i++) {
1460        lower[i]=0.0;
1461        upper[i]=0.0;
1462      }
1463      // ***** Remember to convert to C notation
1464      while (fgets(temp,100,fp)) {
1465        int row;
1466        int value;
1467        if (!strncmp(temp,"ARCS",4))
1468          break;
1469        sscanf(temp,"%d %d",&row,&value);
1470        upper[row-1]=-value;
1471        lower[row-1]=-value;
1472      }
1473      if (strncmp(temp,"ARCS",4)) {
1474        fprintf(stderr,"Unable to find ARCS\n");
1475        exit(2);
1476      }
1477      // number of columns may be underestimate
1478      int * head = new int[2*numberColumns];
1479      int * tail = new int[2*numberColumns];
1480      int * ub = new int[2*numberColumns];
1481      int * cost = new int[2*numberColumns];
1482      // ***** Remember to convert to C notation
1483      numberColumns=0;
1484      while (fgets(temp,100,fp)) {
1485        int iHead;
1486        int iTail;
1487        int iUb;
1488        int iCost;
1489        if (!strncmp(temp,"DEMAND",6))
1490          break;
1491        sscanf(temp,"%d %d %d %d",&iHead,&iTail,&iCost,&iUb);
1492        iHead--;
1493        iTail--;
1494        head[numberColumns]=iHead;
1495        tail[numberColumns]=iTail;
1496        ub[numberColumns]=iUb;
1497        cost[numberColumns]=iCost;
1498        numberColumns++;
1499      }
1500      if (strncmp(temp,"DEMAND",6)) {
1501        fprintf(stderr,"Unable to find DEMAND\n");
1502        exit(2);
1503      }
1504      // ***** Remember to convert to C notation
1505      while (fgets(temp,100,fp)) {
1506        int row;
1507        int value;
1508        if (!strncmp(temp,"END",3))
1509          break;
1510        sscanf(temp,"%d %d",&row,&value);
1511        upper[row-1]=value;
1512        lower[row-1]=value;
1513      }
1514      if (strncmp(temp,"END",3)) {
1515        fprintf(stderr,"Unable to find END\n");
1516        exit(2);
1517      }
1518      printf("Problem %d has %d rows and %d columns\n",problem,
1519             numberRows,numberColumns);
1520      fclose(fp);
1521      ClpSimplex  model;
1522      // now build model
1523     
1524      double * objective =new double[numberColumns];
1525      double * lowerColumn = new double[numberColumns];
1526      double * upperColumn = new double[numberColumns];
1527     
1528      double * element = new double [2*numberColumns];
1529      CoinBigIndex * start = new CoinBigIndex [numberColumns+1];
1530      int * row = new int[2*numberColumns];
1531      start[numberColumns]=2*numberColumns;
1532      for (i=0;i<numberColumns;i++) {
1533        start[i]=2*i;
1534        element[2*i]=-1.0;
1535        element[2*i+1]=1.0;
1536        row[2*i]=head[i];
1537        row[2*i+1]=tail[i];
1538        lowerColumn[i]=0.0;
1539        upperColumn[i]=ub[i];
1540        objective[i]=cost[i];
1541      }
1542      // Create Packed Matrix
1543      CoinPackedMatrix matrix;
1544      int * lengths = NULL;
1545      matrix.assignMatrix(true,numberRows,numberColumns,
1546                          2*numberColumns,element,row,start,lengths);
1547      // load model
1548      model.loadProblem(matrix,
1549                        lowerColumn,upperColumn,objective,
1550                        lower,upper);
1551      model.factorization()->maximumPivots(200+model.numberRows()/100);
1552      model.createStatus();
1553      double time1 = CoinCpuTime();
1554      model.dual();
1555      std::cout<<"Network problem, ClpPackedMatrix took "<<CoinCpuTime()-time1<<" seconds"<<std::endl;
1556      ClpPlusMinusOneMatrix * plusMinus = new ClpPlusMinusOneMatrix(matrix);
1557      assert (plusMinus->getIndices()); // would be zero if not +- one
1558      //ClpPlusMinusOneMatrix *plusminus_matrix;
1559
1560      //plusminus_matrix = new ClpPlusMinusOneMatrix;
1561
1562      //plusminus_matrix->passInCopy(numberRows, numberColumns, true, plusMinus->getMutableIndices(),
1563      //                         plusMinus->startPositive(),plusMinus->startNegative());
1564      model.loadProblem(*plusMinus,
1565                        lowerColumn,upperColumn,objective,
1566                        lower,upper);
1567      //model.replaceMatrix( plusminus_matrix , true);
1568      delete plusMinus;
1569      //model.createStatus();
1570      //model.initialSolve();
1571      //model.writeMps("xx.mps");
1572     
1573      model.factorization()->maximumPivots(200+model.numberRows()/100);
1574      model.createStatus();
1575      time1 = CoinCpuTime();
1576      model.dual();
1577      std::cout<<"Network problem, ClpPlusMinusOneMatrix took "<<CoinCpuTime()-time1<<" seconds"<<std::endl;
1578      ClpNetworkMatrix network(numberColumns,head,tail);
1579      model.loadProblem(network,
1580                        lowerColumn,upperColumn,objective,
1581                        lower,upper);
1582     
1583      model.factorization()->maximumPivots(200+model.numberRows()/100);
1584      model.createStatus();
1585      time1 = CoinCpuTime();
1586      model.dual();
1587      std::cout<<"Network problem, ClpNetworkMatrix took "<<CoinCpuTime()-time1<<" seconds"<<std::endl;
1588      delete [] lower;
1589      delete [] upper;
1590      delete [] head;
1591      delete [] tail;
1592      delete [] ub;
1593      delete [] cost;
1594      delete [] objective;
1595      delete [] lowerColumn;
1596      delete [] upperColumn;
1597    }
1598  }
1599#ifdef QUADRATIC
1600  // Test quadratic to solve linear
1601  if (1) {   
1602    CoinMpsIO m;
1603    std::string fn = mpsDir+"exmip1";
1604    m.readMps(fn.c_str(),"mps");
1605    ClpSimplex solution;
1606    solution.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
1607                         m.getObjCoefficients(),
1608                         m.getRowLower(),m.getRowUpper());
1609    //solution.dual();
1610    // get quadratic part
1611    int numberColumns=solution.numberColumns();
1612    int * start=new int [numberColumns+1];
1613    int * column = new int[numberColumns];
1614    double * element = new double[numberColumns];
1615    int i;
1616    start[0]=0;
1617    int n=0;
1618    int kk=numberColumns-1;
1619    int kk2=numberColumns-1;
1620    for (i=0;i<numberColumns;i++) {
1621      if (i>=kk) {
1622        column[n]=i;
1623        if (i>=kk2)
1624          element[n]=1.0e-1;
1625        else
1626          element[n]=0.0;
1627        n++;
1628      }
1629      start[i+1]=n;
1630    }
1631    // Load up objective
1632    solution.loadQuadraticObjective(numberColumns,start,column,element);
1633    delete [] start;
1634    delete [] column;
1635    delete [] element;
1636    //solution.quadraticSLP(50,1.0e-4);
1637    double objValue = solution.getObjValue();
1638    CoinRelFltEq eq(1.0e-4);
1639    //assert(eq(objValue,820.0));
1640    //solution.setLogLevel(63);
1641    solution.primal();
1642    int numberRows = solution.numberRows();
1643
1644    double * rowPrimal = solution.primalRowSolution();
1645    double * rowDual = solution.dualRowSolution();
1646    double * rowLower = solution.rowLower();
1647    double * rowUpper = solution.rowUpper();
1648   
1649    int iRow;
1650    printf("Rows\n");
1651    for (iRow=0;iRow<numberRows;iRow++) {
1652      printf("%d primal %g dual %g low %g up %g\n",
1653             iRow,rowPrimal[iRow],rowDual[iRow],
1654             rowLower[iRow],rowUpper[iRow]);
1655    }
1656    double * columnPrimal = solution.primalColumnSolution();
1657    double * columnDual = solution.dualColumnSolution();
1658    double * columnLower = solution.columnLower();
1659    double * columnUpper = solution.columnUpper();
1660    objValue = solution.getObjValue();
1661    int iColumn;
1662    printf("Columns\n");
1663    for (iColumn=0;iColumn<numberColumns;iColumn++) {
1664      printf("%d primal %g dual %g low %g up %g\n",
1665             iColumn,columnPrimal[iColumn],columnDual[iColumn],
1666             columnLower[iColumn],columnUpper[iColumn]);
1667    }
1668    //assert(eq(objValue,3.2368421));
1669    //exit(77);
1670  }
1671  // Test quadratic
1672  if (1) {   
1673    CoinMpsIO m;
1674    std::string fn = mpsDir+"share2qp";
1675    //fn = "share2qpb";
1676    m.readMps(fn.c_str(),"mps");
1677    ClpSimplex model;
1678    model.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
1679                         m.getObjCoefficients(),
1680                         m.getRowLower(),m.getRowUpper());
1681    model.dual();
1682    // get quadratic part
1683    int * start=NULL;
1684    int * column = NULL;
1685    double * element = NULL;
1686    m.readQuadraticMps(NULL,start,column,element,2);
1687    int column2[200];
1688    double element2[200];
1689    int start2[80];
1690    int j;
1691    start2[0]=0;
1692    int nel=0;
1693    bool good=false;
1694    for (j=0;j<79;j++) {
1695      if (start[j]==start[j+1]) {
1696        column2[nel]=j;
1697        element2[nel]=0.0;
1698        nel++;
1699      } else {
1700        int i;
1701        for (i=start[j];i<start[j+1];i++) {
1702          column2[nel]=column[i];
1703          element2[nel++]=element[i];
1704        }
1705      }
1706      start2[j+1]=nel;
1707    }
1708    // Load up objective
1709    if (good)
1710      model.loadQuadraticObjective(model.numberColumns(),start2,column2,element2);
1711    else
1712      model.loadQuadraticObjective(model.numberColumns(),start,column,element);
1713    delete [] start;
1714    delete [] column;
1715    delete [] element;
1716    int numberColumns=model.numberColumns();
1717#if 0
1718    model.nonlinearSLP(50,1.0e-4);
1719#else
1720    // Get feasible
1721    ClpObjective * saveObjective = model.objectiveAsObject()->clone();
1722    ClpLinearObjective zeroObjective(NULL,numberColumns);
1723    model.setObjective(&zeroObjective);
1724    model.dual();
1725    model.setObjective(saveObjective);
1726    delete saveObjective;
1727#endif
1728    //model.setLogLevel(63);
1729    //exit(77);
1730    model.setFactorizationFrequency(10);
1731    model.primal();
1732    double objValue = model.getObjValue();
1733    CoinRelFltEq eq(1.0e-4);
1734    assert(eq(objValue,-400.92));
1735  }
1736  if (0) {   
1737    CoinMpsIO m;
1738    std::string fn = "./beale";
1739    //fn = "./jensen";
1740    m.readMps(fn.c_str(),"mps");
1741    ClpSimplex solution;
1742    solution.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
1743                         m.getObjCoefficients(),
1744                         m.getRowLower(),m.getRowUpper());
1745    solution.setDblParam(ClpObjOffset,m.objectiveOffset());
1746    solution.dual();
1747    // get quadratic part
1748    int * start=NULL;
1749    int * column = NULL;
1750    double * element = NULL;
1751    m.readQuadraticMps(NULL,start,column,element,2);
1752    // Load up objective
1753    solution.loadQuadraticObjective(solution.numberColumns(),start,column,element);
1754    delete [] start;
1755    delete [] column;
1756    delete [] element;
1757    solution.primal(1);
1758    solution.nonlinearSLP(50,1.0e-4);
1759    double objValue = solution.getObjValue();
1760    CoinRelFltEq eq(1.0e-4);
1761    assert(eq(objValue,0.5));
1762    solution.primal();
1763    objValue = solution.getObjValue();
1764    assert(eq(objValue,0.5));
1765  }
1766#endif 
1767}
1768void CbcClpUnitTest (const CbcModel & saveModel)
1769{
1770  unsigned int m ;
1771  // See if files exist
1772  FILE * fp;
1773  bool doTest=false;
1774  const char dirsep =  CoinFindDirSeparator();
1775 
1776  // Set directory containing miplib data files.
1777  std::string miplibDir;
1778  miplibDir = dirsep == '/' ? "../../Data/miplib3/" : "..\\..\\Data\\miplib3\\";
1779  std::string test1 = miplibDir +"p0033";
1780  fp=fopen(test1.c_str(),"r");
1781  if (fp) {
1782    doTest=true;
1783    fclose(fp);
1784  }
1785#ifdef COIN_HAS_ZLIB
1786  test1 += ".gz";
1787  fp=fopen(test1.c_str(),"r");
1788  if (fp) {
1789    doTest=true;
1790    fclose(fp);
1791  }
1792#endif
1793  if (!doTest) {
1794    printf("Not doing miplib run as can't find mps files - ? .gz without libz\n");
1795    return;
1796  }
1797  /*
1798    Vectors to hold test problem names and characteristics. The objective value
1799    after optimization (objValue) must agree to the specified tolerance
1800    (objValueTol).
1801  */
1802  std::vector<std::string> mpsName ;
1803  std::vector<int> nRows ;
1804  std::vector<int> nCols ;
1805  std::vector<double> objValueC ;
1806  std::vector<double> objValue ;
1807  std::vector<int> strategy ;
1808  /*
1809    And a macro to make the vector creation marginally readable.
1810  */
1811#define PUSH_MPS(zz_mpsName_zz,\
1812                 zz_nRows_zz,zz_nCols_zz,zz_objValue_zz,zz_objValueC_zz, \
1813                 zz_strategy_zz) \
1814  mpsName.push_back(zz_mpsName_zz) ; \
1815  nRows.push_back(zz_nRows_zz) ; \
1816  nCols.push_back(zz_nCols_zz) ; \
1817  objValueC.push_back(zz_objValueC_zz) ; \
1818  strategy.push_back(zz_strategy_zz) ; \
1819  objValue.push_back(zz_objValue_zz) ;
1820 
1821  /*
1822    Load up the problem vector. Note that the row counts here include the
1823    objective function.
1824   
1825  */
1826  // 0 for no test, 1 for some, 2 for many, 3 for all
1827#define HOWMANY 2
1828#if HOWMANY
1829#if HOWMANY>1
1830  PUSH_MPS("10teams",230,2025,924,917,7);
1831#endif
1832  PUSH_MPS("air03",124,10757,340160,338864.25,7);
1833#if HOWMANY==3
1834  PUSH_MPS("air04",823,8904,56137,55535.436,8);
1835  PUSH_MPS("air05",426,7195,26374,25877.609,8);
1836#endif
1837  //    PUSH_MPS("arki001",1048,1388,7580813.0459,7579599.80787,7);
1838  PUSH_MPS("bell3a",123,133,878430.32,862578.64,7);
1839#if HOWMANY>1
1840  PUSH_MPS("bell5",91,104,8966406.49,8608417.95,7);
1841#endif
1842  PUSH_MPS("blend2",274,353,7.598985,6.9156751140,7);
1843#if HOWMANY>1
1844  PUSH_MPS("cap6000",2176,6000,-2451377,-2451537.325,7);
1845#endif
1846  //    PUSH_MPS("dano3mip",3202,13873,728.1111,576.23162474,7);
1847  //PUSH_MPS("danoint",664,521,65.67,62.637280418,7);
1848  PUSH_MPS("dcmulti",290,548,188182,183975.5397,7);
1849  PUSH_MPS("dsbmip",1182,1886,-305.19817501,-305.19817501,7);
1850  PUSH_MPS("egout",98,141,568.101,149.589,7);
1851  PUSH_MPS("enigma",21,100,0.0,0.0,7);
1852#if HOWMANY==3
1853  PUSH_MPS("fast0507",507,63009,174,172.14556668,7);
1854#endif
1855  PUSH_MPS("fiber",363,1298,405935.18000,156082.51759,7);
1856#if HOWMANY>1
1857  PUSH_MPS("fixnet6",478,878,3983,1200.88,7);
1858#endif
1859  PUSH_MPS("flugpl",18,18,1201500,1167185.7,7);
1860  PUSH_MPS("gen",780,870,112313,112130.0,7);
1861#if HOWMANY>1
1862  PUSH_MPS("gesa2",1392,1224,25779856.372,25476489.678,7);
1863  PUSH_MPS("gesa2_o",1248,1224,25779856.372,25476489.678,7);
1864#endif
1865  PUSH_MPS("gesa3",1368,1152,27991042.648,27833632.451,7);
1866  PUSH_MPS("gesa3_o",1224,1152,27991042.648,27833632.451,7);
1867  PUSH_MPS("gt2",29,188,21166.000,13460.233074,7);
1868#if HOWMANY==3
1869  PUSH_MPS("harp2",112,2993,-73899798.00,-74353341.502,7);
1870#endif
1871  PUSH_MPS("khb05250",101,1350,106940226,95919464.0,7);
1872#if HOWMANY>1
1873  PUSH_MPS("l152lav",97,1989,4722,4656.36,7);
1874#endif
1875  PUSH_MPS("lseu",28,89,1120,834.68,7);
1876  PUSH_MPS("misc03",96,160,3360,1910.,7);
1877  PUSH_MPS("misc06",820,1808,12850.8607,12841.6,7);
1878#if HOWMANY>1
1879  PUSH_MPS("misc07",212,260,2810,1415.0,7);
1880  PUSH_MPS("mitre",2054,10724,115155,114740.5184,7);
1881#endif
1882  PUSH_MPS("mod008",6,319,307,290.9,7);
1883  PUSH_MPS("mod010",146,2655,6548,6532.08,7);
1884#if HOWMANY==3
1885  PUSH_MPS("mod011",4480,10958,-54558535,-62121982.55,7);
1886  PUSH_MPS("modglob",291,422,20740508,20430947.,7);
1887  PUSH_MPS("noswot",182,128,-43,-43.0,7);
1888#endif
1889#if HOWMANY>1
1890  PUSH_MPS("nw04",36,87482,16862,16310.66667,7);
1891#endif
1892  PUSH_MPS("p0033",16,33,3089,2520.57,7);
1893  PUSH_MPS("p0201",133,201,7615,6875.0,7);
1894  PUSH_MPS("p0282",241,282,258411,176867.50,7);
1895  PUSH_MPS("p0548",176,548,8691,315.29,7);
1896  PUSH_MPS("p2756",755,2756,3124,2688.75,7);
1897#if HOWMANY==3
1898  PUSH_MPS("pk1",45,86,11.0,0.0,7);
1899#endif
1900#if HOWMANY>1
1901  PUSH_MPS("pp08a",136,240,7350.0,2748.3452381,7);
1902  PUSH_MPS("pp08aCUTS",246,240,7350.0,5480.6061563,7);
1903#endif
1904#if HOWMANY==3
1905  PUSH_MPS("qiu",1192,840,-132.873137,-931.638857,7);
1906#endif
1907  PUSH_MPS("qnet1",503,1541,16029.692681,14274.102667,7);
1908  PUSH_MPS("qnet1_o",456,1541,16029.692681,12095.571667,7);
1909  PUSH_MPS("rentacar",6803,9557,30356761,28806137.644,7);
1910  PUSH_MPS("rgn",24,180,82.1999,48.7999,7);
1911#if HOWMANY==3
1912  PUSH_MPS("rout",291,556,1077.56,981.86428571,7);
1913  PUSH_MPS("set1ch",492,712,54537.75,32007.73,7);
1914#endif
1915  //    PUSH_MPS("seymour",4944,1372,423,403.84647413,7);
1916  PUSH_MPS("stein27",118,27,18,13.0,7);
1917#if HOWMANY>1
1918  PUSH_MPS("stein45",331,45,30,22.0,7);
1919#endif
1920  PUSH_MPS("vpm1",234,378,20,15.4167,7);
1921  PUSH_MPS("vpm2",234,378,13.75,9.8892645972,7);
1922#endif
1923#undef PUSH_MPS
1924   
1925  int numProbSolved = 0;
1926  double timeTaken=0.0;
1927 
1928  /*
1929    Open the main loop to step through the MPS problems.
1930  */
1931  for (m = 0 ; m < mpsName.size() ; m++) {
1932    std::cerr << "  processing mps file: " << mpsName[m] 
1933              << " (" << m+1 << " out of " << mpsName.size() << ")" << std::endl ;
1934    /*
1935      Stage 1: Read the MPS
1936      and make sure the size of the constraint matrix is correct.
1937    */
1938    CbcModel * model = new CbcModel(saveModel);
1939     
1940    std::string fn = miplibDir+mpsName[m] ;
1941    model->solver()->readMps(fn.c_str(),"") ;
1942    int nr = model->getNumRows() ;
1943    int nc = model->getNumCols() ;
1944    assert(nr == nRows[m]) ;
1945    assert(nc == nCols[m]) ;
1946    /*
1947      Stage 2: Call solver to solve the problem.
1948     
1949      then check the return code and objective.
1950     
1951    */
1952
1953    double startTime = CoinCpuTime();
1954    model->setMaximumNodes(200000);
1955    OsiClpSolverInterface * si =
1956      dynamic_cast<OsiClpSolverInterface *>(model->solver()) ;
1957    assert (si != NULL);
1958    // get clp itself
1959    ClpSimplex * modelC = si->getModelPtr();
1960    modelC->tightenPrimalBounds();
1961    model->initialSolve();
1962    if (modelC->dualBound()==1.0e10) {
1963      // user did not set - so modify
1964      // get largest scaled away from bound
1965      double largest=1.0e-12;
1966      int numberRows = modelC->numberRows();
1967      const double * rowPrimal = modelC->primalRowSolution();
1968      const double * rowLower = modelC->rowLower();
1969      const double * rowUpper = modelC->rowUpper();
1970      const double * rowScale = modelC->rowScale();
1971      int iRow;
1972      for (iRow=0;iRow<numberRows;iRow++) {
1973        double value = rowPrimal[iRow];
1974        double above = value-rowLower[iRow];
1975        double below = rowUpper[iRow]-value;
1976        if (rowScale) {
1977          double multiplier = rowScale[iRow];
1978          above *= multiplier;
1979          below *= multiplier;
1980        }
1981        if (above<1.0e12)
1982          largest = CoinMax(largest,above);
1983        if (below<1.0e12)
1984          largest = CoinMax(largest,below);
1985      }
1986     
1987      int numberColumns = modelC->numberColumns();
1988      const double * columnPrimal = modelC->primalColumnSolution();
1989      const double * columnLower = modelC->columnLower();
1990      const double * columnUpper = modelC->columnUpper();
1991      const double * columnScale = modelC->columnScale();
1992      int iColumn;
1993      for (iColumn=0;iColumn<numberColumns;iColumn++) {
1994        double value = columnPrimal[iColumn];
1995        double above = value-columnLower[iColumn];
1996        double below = columnUpper[iColumn]-value;
1997        if (columnScale) {
1998          double multiplier = 1.0/columnScale[iColumn];
1999          above *= multiplier;
2000          below *= multiplier;
2001        }
2002        if (above<1.0e12)
2003          largest = CoinMax(largest,above);
2004        if (below<1.0e12)
2005          largest = CoinMax(largest,below);
2006      }
2007      //std::cout<<"Largest (scaled) away from bound "<<largest<<std::endl;
2008      modelC->setDualBound(CoinMax(1.0001e8,CoinMin(1000.0*largest,1.00001e10)));
2009    }
2010    model->setMinimumDrop(min(5.0e-2,
2011                                 fabs(model->getMinimizationObjValue())*1.0e-3+1.0e-4));
2012    if (model->getNumCols()<500)
2013      model->setMaximumCutPassesAtRoot(-100); // always do 100 if possible
2014    else if (model->getNumCols()<5000)
2015      model->setMaximumCutPassesAtRoot(100); // use minimum drop
2016    else
2017      model->setMaximumCutPassesAtRoot(20);
2018    model->branchAndBound();
2019     
2020    double timeOfSolution = CoinCpuTime()-startTime;
2021    // Print more statistics
2022    std::cout<<"Cuts at root node changed objective from "<<model->getContinuousObjective()
2023             <<" to "<<model->rootObjectiveAfterCuts()<<std::endl;
2024    int numberGenerators = model->numberCutGenerators();
2025    for (int iGenerator=0;iGenerator<numberGenerators;iGenerator++) {
2026      CbcCutGenerator * generator = model->cutGenerator(iGenerator);
2027      std::cout<<generator->cutGeneratorName()<<" was tried "
2028               <<generator->numberTimesEntered()<<" times and created "
2029               <<generator->numberCutsInTotal()<<" cuts of which "
2030               <<generator->numberCutsActive()<<" were active after adding rounds of cuts";
2031      if (generator->timing())
2032        std::cout<<" ( "<<generator->timeInCutGenerator()<<" seconds)"<<std::endl;
2033      else
2034        std::cout<<std::endl;
2035    }
2036    if (!model->status()) { 
2037      double soln = model->getObjValue();       
2038      CoinRelFltEq eq(1.0e-3) ;
2039      if (eq(soln,objValue[m])) { 
2040        std::cerr
2041          <<"cbc_clp"<<" "
2042          << soln << " = " << objValue[m] << " ; okay";
2043        numProbSolved++;
2044      } else  { 
2045        std::cerr <<"cbc_clp" <<" " <<soln << " != " <<objValue[m] << "; error=" ;
2046        std::cerr <<fabs(objValue[m] - soln); 
2047      }
2048    } else {
2049      std::cerr << "error; too many nodes" ;
2050    }
2051    std::cerr<<" - took " <<timeOfSolution<<" seconds."<<std::endl; 
2052    timeTaken += timeOfSolution;
2053    delete model;
2054  }
2055  std::cerr
2056    <<"cbc_clp" 
2057    <<" solved " 
2058    <<numProbSolved
2059    <<" out of "
2060    <<objValue.size()
2061    <<" and took "
2062    <<timeTaken
2063    <<" seconds."
2064    <<std::endl;
2065}
Note: See TracBrowser for help on using the repository browser.