source: trunk/Cbc/src/unitTestClp.cpp @ 354

Last change on this file since 354 was 354, checked in by ladanyi, 13 years ago

finished switch to subversion

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