source: branches/BSP/trunk/Clp/src/unitTest.cpp @ 1075

Last change on this file since 1075 was 1075, checked in by ladanyi, 12 years ago

Got rid of dependency on Data/Netlib?

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 69.6 KB
Line 
1// Copyright (C) 2002, International Business Machines
2// Corporation and others.  All Rights Reserved.
3
4#include "ClpConfig.h"
5#include "CoinPragma.hpp"
6#include <cassert>
7#include <cstdio>
8#include <cmath>
9#include <cfloat>
10#include <string>
11#include <iostream>
12
13#include "CoinMpsIO.hpp"
14#include "CoinPackedMatrix.hpp"
15#include "CoinPackedVector.hpp"
16#include "CoinHelperFunctions.hpp"
17#include "CoinTime.hpp"
18
19#include "ClpFactorization.hpp"
20#include "ClpSimplex.hpp"
21#include "ClpSimplexOther.hpp"
22#include "ClpSimplexNonlinear.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}
189static void printSol(ClpSimplex & model) 
190{
191  int numberRows = model.numberRows();
192  int numberColumns = model.numberColumns();
193 
194  double * rowPrimal = model.primalRowSolution();
195  double * rowDual = model.dualRowSolution();
196  double * rowLower = model.rowLower();
197  double * rowUpper = model.rowUpper();
198 
199  int iRow;
200  double objValue = model.getObjValue();
201  printf("Objvalue %g Rows (%d)\n",objValue,numberRows);
202  for (iRow=0;iRow<numberRows;iRow++) {
203    printf("%d primal %g dual %g low %g up %g\n",
204           iRow,rowPrimal[iRow],rowDual[iRow],
205           rowLower[iRow],rowUpper[iRow]);
206  }
207  double * columnPrimal = model.primalColumnSolution();
208  double * columnDual = model.dualColumnSolution();
209  double * columnLower = model.columnLower();
210  double * columnUpper = model.columnUpper();
211  double offset;
212  //const double * gradient = model.objectiveAsObject()->gradient(&model,
213  //                                                       columnPrimal,offset,true,1);
214  const double * gradient = model.objective(columnPrimal,offset);
215  int iColumn;
216  objValue = -offset-model.objectiveOffset();
217  printf("offset %g (%g)\n",offset,model.objectiveOffset());
218  printf("Columns (%d)\n",numberColumns);
219  for (iColumn=0;iColumn<numberColumns;iColumn++) {
220    printf("%d primal %g dual %g low %g up %g\n",
221           iColumn,columnPrimal[iColumn],columnDual[iColumn],
222           columnLower[iColumn],columnUpper[iColumn]);
223    objValue += columnPrimal[iColumn]*gradient[iColumn];
224    if (fabs(columnPrimal[iColumn]*gradient[iColumn])>1.0e-8)
225      printf("obj -> %g gradient %g\n",objValue,gradient[iColumn]);
226  }
227  printf("Computed objective %g\n",objValue);
228}
229//----------------------------------------------------------------
230// unitTest [-mpsDir=V1] [-netlibDir=V2] [-test]
231//
232// where:
233//   -mpsDir: directory containing mps test files
234//       Default value V1="../../Data/Sample"   
235//   -netlibDir: directory containing netlib files
236//       Default value V2="../../Data/Netlib"
237//   -test
238//       If specified, then netlib test set run
239//
240// All parameters are optional.
241//----------------------------------------------------------------
242int mainTest (int argc, const char *argv[],int algorithm,
243              ClpSimplex empty, bool doPresolve, int switchOffValue,bool doVector)
244{
245  int i;
246
247  if (switchOffValue>0) {
248    // switch off some
249    int iTest;
250    for (iTest=0;iTest<NUMBER_ALGORITHMS;iTest++) {
251      int bottom = switchOffValue%10;
252      assert (bottom==0||bottom==1);
253      switchOffValue/= 10;
254      switchOff[iTest]=0;
255    }
256  }
257
258  // define valid parameter keywords
259  std::set<std::string> definedKeyWords;
260  definedKeyWords.insert("-mpsDir");
261  definedKeyWords.insert("-netlibDir");
262  definedKeyWords.insert("-netlib");
263
264  // Create a map of parameter keys and associated data
265  std::map<std::string,std::string> parms;
266  for ( i=1; i<argc; i++ ) {
267    std::string parm(argv[i]);
268    std::string key,value;
269    std::string::size_type  eqPos = parm.find('=');
270
271    // Does parm contain and '='
272    if ( eqPos==std::string::npos ) {
273      //Parm does not contain '='
274      key = parm;
275    }
276    else {
277      key=parm.substr(0,eqPos);
278      value=parm.substr(eqPos+1);
279    }
280
281    // Is specifed key valid?
282    if ( definedKeyWords.find(key) == definedKeyWords.end() ) {
283      // invalid key word.
284      // Write help text
285      std::cerr <<"Undefined parameter \"" <<key <<"\".\n";
286      std::cerr <<"Correct usage: \n";
287      std::cerr <<"  unitTest [-mpsDir=V1] [-netlibDir=V2] [-test[=V3]]\n";
288      std::cerr <<"  where:\n";
289      std::cerr <<"    -mpsDir: directory containing mps test files\n";
290      std::cerr <<"        Default value V1=\"../../Data/Sample\"\n";
291      std::cerr <<"    -netlibDir: directory containing netlib files\n";
292      std::cerr <<"        Default value V2=\"../../Data/Netlib\"\n";
293      std::cerr <<"    -test\n";
294      std::cerr <<"        If specified, then netlib testset run.\n";
295      std::cerr <<"        If V3 then taken as single file\n";
296      return 1;
297    }
298    parms[key]=value;
299  }
300 
301  const char dirsep =  CoinFindDirSeparator();
302  // Set directory containing mps data files.
303  std::string mpsDir;
304  if (parms.find("-mpsDir") != parms.end())
305    mpsDir=parms["-mpsDir"] + dirsep;
306  else 
307    mpsDir = dirsep == '/' ? "../../Data/Sample/" : "..\\..\\Data\\Sample\\";
308 
309  // Set directory containing netlib data files.
310  std::string netlibDir;
311  if (parms.find("-netlibDir") != parms.end())
312    netlibDir=parms["-netlibDir"] + dirsep;
313  else 
314    netlibDir = dirsep == '/' ? "../../Data/Netlib/" : "..\\..\\Data\\Netlib\\";
315  if (!empty.numberRows()) {
316    testingMessage( "Testing ClpSimplex\n" );
317    ClpSimplexUnitTest(mpsDir,netlibDir);
318  }
319  if (parms.find("-netlib") != parms.end()||empty.numberRows())
320  {
321    unsigned int m;
322   
323    // Define test problems:
324    //   mps names,
325    //   maximization or minimization,
326    //   Number of rows and columns in problem, and
327    //   objective function value
328    std::vector<std::string> mpsName;
329    std::vector<bool> min;
330    std::vector<int> nRows;
331    std::vector<int> nCols;
332    std::vector<double> objValue;
333    std::vector<double> objValueTol;
334    // 100 added means no presolve
335    std::vector<int> bestStrategy;
336    if(empty.numberRows()) {
337      std::string alg;
338      for (int iTest=0;iTest<NUMBER_ALGORITHMS;iTest++) {
339        ClpSolve solveOptions=setupForSolve(iTest,alg,0);
340        printf("%d %s ",iTest,alg.c_str());
341        if (switchOff[iTest]) 
342          printf("skipped by user\n");
343        else if(solveOptions.getSolveType()==ClpSolve::notImplemented)
344          printf("skipped as not available\n");
345        else
346          printf("will be tested\n");
347      }
348    }
349    if (!empty.numberRows()) {
350      mpsName.push_back("25fv47");
351      min.push_back(true);
352      nRows.push_back(822);
353      nCols.push_back(1571);
354      objValueTol.push_back(1.E-10);
355      objValue.push_back(5.5018458883E+03);
356      bestStrategy.push_back(0);
357      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);
358      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);
359      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);
360      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);
361      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);
362      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);
363      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);
364      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);
365      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);
366      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);
367      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);
368      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);
369      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);
370      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);
371      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);
372      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);
373      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);
374      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);
375      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);
376      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);
377      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);
378      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);
379      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);
380      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);
381      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);
382      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);
383      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);
384      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);
385      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.
386      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);
387      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);
388      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);
389      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);
390      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);
391      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);
392      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);
393      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);
394      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);
395      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);
396      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);
397      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);
398      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);
399      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);
400      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);
401      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);
402      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);
403      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);
404      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);
405      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);
406      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);
407      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);
408      //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);
409      //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);
410      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);
411      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);
412      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);
413      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);
414      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);
415      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);
416      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);
417      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);
418      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);
419      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);
420      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);
421      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);
422      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);
423      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);
424      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);
425      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);
426      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);
427      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);
428      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);
429      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);
430      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);
431      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);
432      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);
433      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);
434      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);
435      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);
436      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);
437      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);
438      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);
439      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);
440      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);
441      //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);
442      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);
443      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);
444      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);
445      //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);
446      //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);
447      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);
448      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);
449      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);
450      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);
451    } else {
452      // Just testing one
453      mpsName.push_back(empty.problemName());min.push_back(true);nRows.push_back(-1);
454      nCols.push_back(-1);objValueTol.push_back(1.e-10);
455      objValue.push_back(0.0);bestStrategy.push_back(0);
456      int iTest;
457      std::string alg;
458      for (iTest=0;iTest<NUMBER_ALGORITHMS;iTest++) {
459        ClpSolve solveOptions=setupForSolve(iTest,alg,0);
460        printf("%d %s ",iTest,alg.c_str());
461        if (switchOff[iTest]) 
462          printf("skipped by user\n");
463        else if(solveOptions.getSolveType()==ClpSolve::notImplemented)
464          printf("skipped as not available\n");
465        else
466          printf("will be tested\n");
467      }
468    }
469
470    double timeTaken =0.0;
471    if( !barrierAvailable)
472      switchOff[0]=1;
473  // Loop once for each Mps File
474    for (m=0; m<mpsName.size(); m++ ) {
475      std::cerr <<"  processing mps file: " <<mpsName[m] 
476                <<" (" <<m+1 <<" out of " <<mpsName.size() <<")" <<std::endl;
477
478      ClpSimplex solutionBase=empty;
479      std::string fn = netlibDir+mpsName[m];
480      if (!empty.numberRows()||algorithm<6) {
481        // Read data mps file,
482        CoinMpsIO mps;
483        int nerrors=mps.readMps(fn.c_str(),"mps");
484        if (nerrors) {
485          std::cerr << "Error " << nerrors << " when reading model from "
486                    << fn.c_str() << "! Aborting tests.\n";
487          return 1;
488        }
489        solutionBase.loadProblem(*mps.getMatrixByCol(),mps.getColLower(),
490                                 mps.getColUpper(),
491                                 mps.getObjCoefficients(),
492                                 mps.getRowLower(),mps.getRowUpper());
493       
494        solutionBase.setDblParam(ClpObjOffset,mps.objectiveOffset());
495      } 
496     
497      // Runs through strategies
498      if (algorithm==6||algorithm==7) {
499        // algorithms tested are at top of file
500        double testTime[NUMBER_ALGORITHMS];
501        std::string alg[NUMBER_ALGORITHMS];
502        int iTest;
503        for (iTest=0;iTest<NUMBER_ALGORITHMS;iTest++) {
504          ClpSolve solveOptions=setupForSolve(iTest,alg[iTest],1);
505          if (solveOptions.getSolveType()!=ClpSolve::notImplemented) {
506            double time1 = CoinCpuTime();
507            ClpSimplex solution=solutionBase;
508            if (solution.maximumSeconds()<0.0)
509              solution.setMaximumSeconds(120.0);
510            if (doVector) {
511              ClpMatrixBase * matrix = solution.clpMatrix();
512              if (dynamic_cast< ClpPackedMatrix*>(matrix)) {
513                ClpPackedMatrix * clpMatrix = dynamic_cast< ClpPackedMatrix*>(matrix);
514                clpMatrix->makeSpecialColumnCopy();
515              }
516            }
517            solution.initialSolve(solveOptions);
518            double time2 = CoinCpuTime()-time1;
519            testTime[iTest]=time2;
520            printf("Took %g seconds - status %d\n",time2,solution.problemStatus());
521            if (solution.problemStatus()) 
522              testTime[iTest]=1.0e20;
523          } else {
524            testTime[iTest]=1.0e30;
525          }
526        }
527        int iBest=-1;
528        double dBest=1.0e10;
529        printf("%s",fn.c_str());
530        for (iTest=0;iTest<NUMBER_ALGORITHMS;iTest++) {
531          if (testTime[iTest]<1.0e30) {
532            printf(" %s %g",
533                   alg[iTest].c_str(),testTime[iTest]);
534            if (testTime[iTest]<dBest) {
535              dBest=testTime[iTest];
536              iBest=iTest;
537            }
538          }
539        }
540        printf("\n");
541        if (iBest>=0)
542          printf("Best strategy for %s is %s (%d) which takes %g seconds\n",
543                 fn.c_str(),alg[iBest].c_str(),iBest,testTime[iBest]);
544        else
545          printf("No strategy finished in time limit\n");
546        continue;
547      }
548      double time1 = CoinCpuTime();
549      ClpSimplex solution=solutionBase;
550#if 0
551      solution.setOptimizationDirection(-1);
552      {
553        int j;
554        double * obj = solution.objective();
555        int n=solution.numberColumns();
556        for (j=0;j<n;j++)
557          obj[j] *= -1.0;
558      }
559#endif
560      ClpSolve::SolveType method;
561      ClpSolve::PresolveType presolveType;
562      ClpSolve solveOptions;
563      if (doPresolve)
564        presolveType=ClpSolve::presolveOn;
565      else
566        presolveType=ClpSolve::presolveOff;
567      solveOptions.setPresolveType(presolveType,5);
568      std::string nameAlgorithm;
569      if (algorithm!=5) {
570        if (algorithm==0) {
571          method=ClpSolve::useDual;
572          nameAlgorithm="dual";
573        } else if (algorithm==1) {
574          method=ClpSolve::usePrimalorSprint;
575          nameAlgorithm="primal";
576        } else if (algorithm==3) {
577          method=ClpSolve::automatic;
578          nameAlgorithm="either";
579        } else {
580          nameAlgorithm="barrier-slow";
581#ifdef UFL_BARRIER
582          solveOptions.setSpecialOption(4,4);
583          nameAlgorithm="barrier-UFL";
584#endif
585#ifdef WSSMP_BARRIER
586          solveOptions.setSpecialOption(4,2);
587          nameAlgorithm="barrier-WSSMP";
588#endif
589          method = ClpSolve::useBarrier;
590        }
591        solveOptions.setSolveType(method);
592      } else {
593        int iAlg = bestStrategy[m];
594        int presolveOff=iAlg/100;
595        iAlg=iAlg % 100;
596        if( !barrierAvailable&&iAlg==0) {
597          if (nRows[m]!=2172)
598            iAlg = 5; // try primal
599          else
600            iAlg=3; // d2q06c
601        }
602        solveOptions=setupForSolve(iAlg,nameAlgorithm,0);
603        if (presolveOff)
604          solveOptions.setPresolveType(ClpSolve::presolveOff);
605      }
606      if (doVector) {
607        ClpMatrixBase * matrix = solution.clpMatrix();
608        if (dynamic_cast< ClpPackedMatrix*>(matrix)) {
609          ClpPackedMatrix * clpMatrix = dynamic_cast< ClpPackedMatrix*>(matrix);
610          clpMatrix->makeSpecialColumnCopy();
611        }
612      }
613      solution.initialSolve(solveOptions);
614      double time2 = CoinCpuTime()-time1;
615      timeTaken += time2;
616      printf("%s took %g seconds using algorithm %s\n",fn.c_str(),time2,nameAlgorithm.c_str());
617      // Test objective solution value
618      {
619        double soln = solution.objectiveValue();
620        CoinRelFltEq eq(objValueTol[m]);
621        std::cerr <<soln <<",  " <<objValue[m] <<" diff "<<
622          soln-objValue[m]<<std::endl;
623        if(!eq(soln,objValue[m]))
624          printf("** difference fails\n");
625      }
626    }
627    printf("Total time %g seconds\n",timeTaken);
628  }
629  else {
630    testingMessage( "***Skipped Testing on netlib    ***\n" );
631    testingMessage( "***use -netlib to test class***\n" );
632  }
633 
634  testingMessage( "All tests completed successfully\n" );
635  return 0;
636}
637
638 
639// Display message on stdout and stderr
640void testingMessage( const char * const msg )
641{
642  std::cerr <<msg;
643  //cout <<endl <<"*****************************************"
644  //     <<endl <<msg <<endl;
645}
646
647//--------------------------------------------------------------------------
648// test factorization methods and simplex method and simple barrier
649void
650ClpSimplexUnitTest(const std::string & mpsDir,
651                   const std::string & netlibDir)
652{
653 
654  CoinRelFltEq eq(0.000001);
655
656  {
657    ClpSimplex solution;
658 
659    // matrix data
660    //deliberate hiccup of 2 between 0 and 1
661    CoinBigIndex start[5]={0,4,7,8,9};
662    int length[5]={2,3,1,1,1};
663    int rows[11]={0,2,-1,-1,0,1,2,0,1,2};
664    double elements[11]={7.0,2.0,1.0e10,1.0e10,-2.0,1.0,-2.0,1,1,1};
665    CoinPackedMatrix matrix(true,3,5,8,elements,rows,start,length);
666   
667    // rim data
668    double objective[7]={-4.0,1.0,0.0,0.0,0.0,0.0,0.0};
669    double rowLower[5]={14.0,3.0,3.0,1.0e10,1.0e10};
670    double rowUpper[5]={14.0,3.0,3.0,-1.0e10,-1.0e10};
671    double colLower[7]={0.0,0.0,0.0,0.0,0.0,0.0,0.0};
672    double colUpper[7]={100.0,100.0,100.0,100.0,100.0,100.0,100.0};
673   
674    // basis 1
675    int rowBasis1[3]={-1,-1,-1};
676    int colBasis1[5]={1,1,-1,-1,1};
677    solution.loadProblem(matrix,colLower,colUpper,objective,
678                         rowLower,rowUpper);
679    int i;
680    solution.createStatus();
681    for (i=0;i<3;i++) {
682      if (rowBasis1[i]<0) {
683        solution.setRowStatus(i,ClpSimplex::atLowerBound);
684      } else {
685        solution.setRowStatus(i,ClpSimplex::basic);
686      }
687    }
688    for (i=0;i<5;i++) {
689      if (colBasis1[i]<0) {
690        solution.setColumnStatus(i,ClpSimplex::atLowerBound);
691      } else {
692        solution.setColumnStatus(i,ClpSimplex::basic);
693      }
694    }
695    solution.setLogLevel(3+4+8+16+32);
696    solution.primal();
697    for (i=0;i<3;i++) {
698      if (rowBasis1[i]<0) {
699        solution.setRowStatus(i,ClpSimplex::atLowerBound);
700      } else {
701        solution.setRowStatus(i,ClpSimplex::basic);
702      }
703    }
704    for (i=0;i<5;i++) {
705      if (colBasis1[i]<0) {
706        solution.setColumnStatus(i,ClpSimplex::atLowerBound);
707      } else {
708        solution.setColumnStatus(i,ClpSimplex::basic);
709      }
710    }
711    // intricate stuff does not work with scaling
712    solution.scaling(0);
713    int returnCode = solution.factorize ( );
714    assert(!returnCode);
715    const double * colsol = solution.primalColumnSolution();
716    const double * rowsol = solution.primalRowSolution();
717    solution.getSolution(rowsol,colsol);
718    double colsol1[5]={20.0/7.0,3.0,0.0,0.0,23.0/7.0};
719    for (i=0;i<5;i++) {
720      assert(eq(colsol[i],colsol1[i]));
721    }
722    // now feed in again without actually doing factorization
723    ClpFactorization factorization2 = *solution.factorization();
724    ClpSimplex solution2 = solution;
725    solution2.setFactorization(factorization2);
726    solution2.createStatus();
727    for (i=0;i<3;i++) {
728      if (rowBasis1[i]<0) {
729        solution2.setRowStatus(i,ClpSimplex::atLowerBound);
730      } else {
731        solution2.setRowStatus(i,ClpSimplex::basic);
732      }
733    }
734    for (i=0;i<5;i++) {
735      if (colBasis1[i]<0) {
736        solution2.setColumnStatus(i,ClpSimplex::atLowerBound);
737      } else {
738        solution2.setColumnStatus(i,ClpSimplex::basic);
739      }
740    }
741    // intricate stuff does not work with scaling
742    solution2.scaling(0);
743    solution2.getSolution(rowsol,colsol);
744    colsol = solution2.primalColumnSolution();
745    rowsol = solution2.primalRowSolution();
746    for (i=0;i<5;i++) {
747      assert(eq(colsol[i],colsol1[i]));
748    }
749    solution2.setDualBound(0.1);
750    solution2.dual();
751    objective[2]=-1.0;
752    objective[3]=-0.5;
753    objective[4]=10.0;
754    solution.dual();
755    for (i=0;i<3;i++) {
756      rowLower[i]=-1.0e20;
757      colUpper[i+2]=0.0;
758    }
759    solution.setLogLevel(3);
760    solution.dual();
761    double rowObjective[]={1.0,0.5,-10.0};
762    solution.loadProblem(matrix,colLower,colUpper,objective,
763                         rowLower,rowUpper,rowObjective);
764    solution.dual();
765    solution.loadProblem(matrix,colLower,colUpper,objective,
766                         rowLower,rowUpper,rowObjective);
767    solution.primal();
768  }
769#ifndef COIN_NO_CLP_MESSAGE
770  {   
771    CoinMpsIO m;
772    std::string fn = mpsDir+"exmip1";
773    m.readMps(fn.c_str(),"mps");
774    ClpSimplex solution;
775    solution.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
776                         m.getObjCoefficients(),
777                         m.getRowLower(),m.getRowUpper());
778    solution.dual();
779    // Test event handling
780    MyEventHandler handler;
781    solution.passInEventHandler(&handler);
782    int numberRows=solution.numberRows();
783    // make sure values pass has something to do
784    for (int i=0;i<numberRows;i++)
785      solution.setRowStatus(i,ClpSimplex::basic);
786    solution.primal(1);
787    assert (solution.secondaryStatus()==102); // Came out at end of pass
788  }
789  // Test Message handler
790  {   
791    CoinMpsIO m;
792    std::string fn = mpsDir+"exmip1";
793    //fn = "Test/subGams4";
794    m.readMps(fn.c_str(),"mps");
795    ClpSimplex model;
796    model.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
797                         m.getObjCoefficients(),
798                         m.getRowLower(),m.getRowUpper());
799    // Message handler
800    MyMessageHandler messageHandler(&model);
801    std::cout<<"Testing derived message handler"<<std::endl;
802    model.passInMessageHandler(&messageHandler);
803    model.messagesPointer()->setDetailMessage(1,102);
804    model.setFactorizationFrequency(10);
805    model.primal();
806    model.primal(0,3);
807    model.setObjCoeff(3,-2.9473684210526314);
808    model.primal(0,3);
809    // Write saved solutions
810    int nc = model.getNumCols();
811    int s; 
812    std::deque<StdVectorDouble> fep = messageHandler.getFeasibleExtremePoints();
813    int numSavedSolutions = fep.size();
814    for ( s=0; s<numSavedSolutions; ++s ) {
815      const StdVectorDouble & solnVec = fep[s];
816      for ( int c=0; c<nc; ++c ) {
817        if (fabs(solnVec[c])>1.0e-8)
818          std::cout <<"Saved Solution: " <<s <<" ColNum: " <<c <<" Value: " <<solnVec[c] <<std::endl;
819      }
820    }
821    // Solve again without scaling
822    // and maximize then minimize
823    messageHandler.clearFeasibleExtremePoints();
824    model.scaling(0);
825    model.setOptimizationDirection(-1);
826    model.primal();
827    model.setOptimizationDirection(1);
828    model.primal();
829    fep = messageHandler.getFeasibleExtremePoints();
830    numSavedSolutions = fep.size();
831    for ( s=0; s<numSavedSolutions; ++s ) {
832      const StdVectorDouble & solnVec = fep[s];
833      for ( int c=0; c<nc; ++c ) {
834        if (fabs(solnVec[c])>1.0e-8)
835          std::cout <<"Saved Solution: " <<s <<" ColNum: " <<c <<" Value: " <<solnVec[c] <<std::endl;
836      }
837    }
838  }
839#endif
840  // Test dual ranging
841  {   
842    CoinMpsIO m;
843    std::string fn = mpsDir+"exmip1";
844    m.readMps(fn.c_str(),"mps");
845    ClpSimplex model;
846    model.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
847                         m.getObjCoefficients(),
848                         m.getRowLower(),m.getRowUpper());
849    model.primal();
850    int which[13] = {0,1,2,3,4,5,6,7,8,9,10,11,12};
851    double costIncrease[13];
852    int sequenceIncrease[13];
853    double costDecrease[13];
854    int sequenceDecrease[13];
855    // ranging
856    model.dualRanging(13,which,costIncrease,sequenceIncrease,
857                      costDecrease,sequenceDecrease);
858    int i;
859    for ( i=0;i<13;i++)
860      printf("%d increase %g %d, decrease %g %d\n",
861             i,costIncrease[i],sequenceIncrease[i],
862             costDecrease[i],sequenceDecrease[i]);
863    assert (fabs(costDecrease[3])<1.0e-4);
864    assert (fabs(costIncrease[7]-1.0)<1.0e-4);
865    model.setOptimizationDirection(-1);
866    {
867      int j;
868      double * obj = model.objective();
869      int n=model.numberColumns();
870      for (j=0;j<n;j++) 
871        obj[j] *= -1.0;
872    }
873    double costIncrease2[13];
874    int sequenceIncrease2[13];
875    double costDecrease2[13];
876    int sequenceDecrease2[13];
877    // ranging
878    model.dualRanging(13,which,costIncrease2,sequenceIncrease2,
879                      costDecrease2,sequenceDecrease2);
880    for (i=0;i<13;i++) {
881      assert (fabs(costIncrease[i]-costDecrease2[i])<1.0e-6);
882      assert (fabs(costDecrease[i]-costIncrease2[i])<1.0e-6);
883      assert (sequenceIncrease[i]==sequenceDecrease2[i]);
884      assert (sequenceDecrease[i]==sequenceIncrease2[i]);
885    }
886    // Now delete all rows and see what happens
887    model.deleteRows(model.numberRows(),which);
888    model.primal();
889    // ranging
890    if (!model.dualRanging(8,which,costIncrease,sequenceIncrease,
891                           costDecrease,sequenceDecrease)) {
892      for (i=0;i<8;i++) {
893        printf("%d increase %g %d, decrease %g %d\n",
894               i,costIncrease[i],sequenceIncrease[i],
895               costDecrease[i],sequenceDecrease[i]);
896      }
897    }
898  }
899  // Test primal ranging
900  {   
901    CoinMpsIO m;
902    std::string fn = mpsDir+"exmip1";
903    m.readMps(fn.c_str(),"mps");
904    ClpSimplex model;
905    model.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
906                         m.getObjCoefficients(),
907                         m.getRowLower(),m.getRowUpper());
908    model.primal();
909    int which[13] = {0,1,2,3,4,5,6,7,8,9,10,11,12};
910    double valueIncrease[13];
911    int sequenceIncrease[13];
912    double valueDecrease[13];
913    int sequenceDecrease[13];
914    // ranging
915    model.primalRanging(13,which,valueIncrease,sequenceIncrease,
916                      valueDecrease,sequenceDecrease);
917    int i;
918    for ( i=0;i<13;i++)
919      printf("%d increase %g %d, decrease %g %d\n",
920             i,valueIncrease[i],sequenceIncrease[i],
921             valueDecrease[i],sequenceDecrease[i]);
922    assert (fabs(valueIncrease[3]-0.642857)<1.0e-4);
923    assert (fabs(valueIncrease[8]-2.95113)<1.0e-4);
924#if 0
925    // out until I find optimization bug
926    // Test parametrics
927    ClpSimplexOther * model2 = (ClpSimplexOther *) (&model);
928    double rhs[]={ 1.0,2.0,3.0,4.0,5.0};
929    double endingTheta=1.0;
930    model2->scaling(0);
931    model2->setLogLevel(63);
932    model2->parametrics(0.0,endingTheta,0.1,
933                        NULL,NULL,rhs,rhs,NULL);
934#endif
935  }
936  // Test binv etc
937  {   
938    /*
939       Wolsey : Page 130
940       max 4x1 -  x2
941       7x1 - 2x2    <= 14
942       x2    <= 3
943       2x1 - 2x2    <= 3
944       x1 in Z+, x2 >= 0
945
946       note slacks are -1 in Clp so signs may be different
947    */
948   
949    int n_cols = 2;
950    int n_rows = 3;
951   
952    double obj[2] = {-4.0, 1.0};
953    double collb[2] = {0.0, 0.0};
954    double colub[2] = {COIN_DBL_MAX, COIN_DBL_MAX};
955    double rowlb[3] = {-COIN_DBL_MAX, -COIN_DBL_MAX, -COIN_DBL_MAX};
956    double rowub[3] = {14.0, 3.0, 3.0};
957   
958    int rowIndices[5] =  {0,     2,    0,    1,    2};
959    int colIndices[5] =  {0,     0,    1,    1,    1};
960    double elements[5] = {7.0, 2.0, -2.0,  1.0, -2.0};
961    CoinPackedMatrix M(true, rowIndices, colIndices, elements, 5);
962
963    ClpSimplex model;
964    model.loadProblem(M, collb, colub, obj, rowlb, rowub);
965    model.dual(0,1); // keep factorization
966   
967    //check that the tableau matches wolsey (B-1 A)
968    // slacks in second part of binvA
969    double * binvA = (double*) malloc((n_cols+n_rows) * sizeof(double));
970   
971    printf("B-1 A by row\n");
972    int i;
973    for( i = 0; i < n_rows; i++){
974      model.getBInvARow(i, binvA,binvA+n_cols);
975      printf("row: %d -> ",i);
976      for(int j=0; j < n_cols+n_rows; j++){
977        printf("%g, ", binvA[j]);
978      }
979      printf("\n");
980    }
981    // See if can re-use factorization AND arrays
982    model.primal(0,3+4); // keep factorization
983    // And do by column
984    printf("B-1 A by column\n");
985    for( i = 0; i < n_rows+n_cols; i++){
986      model.getBInvACol(i, binvA);
987      printf("column: %d -> ",i);
988      for(int j=0; j < n_rows; j++){
989        printf("%g, ", binvA[j]);
990      }
991      printf("\n");
992    }
993    free(binvA);
994    model.setColUpper(1,2.0);
995    model.dual(0,2+4); // use factorization and arrays
996    model.dual(0,2); // hopefully will not use factorization
997    model.primal(0,3+4); // keep factorization
998    // but say basis has changed
999    model.setWhatsChanged(model.whatsChanged()&(~512));
1000    model.dual(0,2); // hopefully will not use factorization
1001  }
1002  // test steepest edge
1003  {   
1004    CoinMpsIO m;
1005    std::string fn = mpsDir+"finnis";
1006    int returnCode = m.readMps(fn.c_str(),"mps");
1007    if (returnCode) {
1008      // probable cause is that gz not there
1009      fprintf(stderr,"Unable to open finnis.mps in %s\n", mpsDir.c_str());
1010      fprintf(stderr,"Most probable cause is finnis.mps is gzipped i.e. finnis.mps.gz and libz has not been activated\n");
1011      fprintf(stderr,"Either gunzip files or edit Makefiles/Makefile.location to get libz\n");
1012      exit(999);
1013    }
1014    ClpModel model;
1015    model.loadProblem(*m.getMatrixByCol(),m.getColLower(),
1016                    m.getColUpper(),
1017                    m.getObjCoefficients(),
1018                    m.getRowLower(),m.getRowUpper());
1019    ClpSimplex solution(model);
1020
1021    solution.scaling(1); 
1022    solution.setDualBound(1.0e8);
1023    //solution.factorization()->maximumPivots(1);
1024    //solution.setLogLevel(3);
1025    solution.setDualTolerance(1.0e-7);
1026    // set objective sense,
1027    ClpDualRowSteepest steep;
1028    solution.setDualRowPivotAlgorithm(steep);
1029    solution.setDblParam(ClpObjOffset,m.objectiveOffset());
1030    solution.dual();
1031  }
1032  // test normal solution
1033  {   
1034    CoinMpsIO m;
1035    std::string fn = mpsDir+"afiro";
1036    m.readMps(fn.c_str(),"mps");
1037    ClpSimplex solution;
1038    ClpModel model;
1039    // do twice - without and with scaling
1040    int iPass;
1041    for (iPass=0;iPass<2;iPass++) {
1042      // explicit row objective for testing
1043      int nr = m.getNumRows();
1044      double * rowObj = new double[nr];
1045      CoinFillN(rowObj,nr,0.0);
1046      model.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
1047                      m.getObjCoefficients(),
1048                      m.getRowLower(),m.getRowUpper(),rowObj);
1049      delete [] rowObj;
1050      solution = ClpSimplex(model);
1051      if (iPass) {
1052        solution.scaling();
1053      }
1054      solution.dual();
1055      solution.dual();
1056      // test optimal
1057      assert (solution.status()==0);
1058      int numberColumns = solution.numberColumns();
1059      int numberRows = solution.numberRows();
1060      CoinPackedVector colsol(numberColumns,solution.primalColumnSolution());
1061      double * objective = solution.objective();
1062      double objValue = colsol.dotProduct(objective);
1063      CoinRelFltEq eq(1.0e-8);
1064      assert(eq(objValue,-4.6475314286e+02));
1065      // Test auxiliary model
1066      //solution.scaling(0);
1067      solution.auxiliaryModel(63-2); // bounds may change
1068      solution.dual();
1069      solution.primal();
1070      solution.allSlackBasis();
1071      solution.dual();
1072      assert(eq(solution.objectiveValue(),-4.6475314286e+02));
1073      solution.auxiliaryModel(-1);
1074      solution.dual();
1075      assert(eq(solution.objectiveValue(),-4.6475314286e+02));
1076      double * lower = solution.columnLower();
1077      double * upper = solution.columnUpper();
1078      double * sol = solution.primalColumnSolution();
1079      double * result = new double[numberColumns];
1080      CoinFillN ( result, numberColumns,0.0);
1081      solution.matrix()->transposeTimes(solution.dualRowSolution(), result);
1082      int iRow , iColumn;
1083      // see if feasible and dual feasible
1084      for (iColumn=0;iColumn<numberColumns;iColumn++) {
1085        double value = sol[iColumn];
1086        assert(value<upper[iColumn]+1.0e-8);
1087        assert(value>lower[iColumn]-1.0e-8);
1088        value = objective[iColumn]-result[iColumn];
1089        assert (value>-1.0e-5);
1090        if (sol[iColumn]>1.0e-5)
1091          assert (value<1.0e-5);
1092      }
1093      delete [] result;
1094      result = new double[numberRows];
1095      CoinFillN ( result, numberRows,0.0);
1096      solution.matrix()->times(colsol, result);
1097      lower = solution.rowLower();
1098      upper = solution.rowUpper();
1099      sol = solution.primalRowSolution();
1100      for (iRow=0;iRow<numberRows;iRow++) {
1101        double value = result[iRow];
1102        assert(eq(value,sol[iRow]));
1103        assert(value<upper[iRow]+1.0e-8);
1104        assert(value>lower[iRow]-1.0e-8);
1105      }
1106      delete [] result;
1107      // test row objective
1108      double * rowObjective = solution.rowObjective();
1109      CoinDisjointCopyN(solution.dualRowSolution(),numberRows,rowObjective);
1110      CoinDisjointCopyN(solution.dualColumnSolution(),numberColumns,objective);
1111      // this sets up all slack basis
1112      solution.createStatus();
1113      solution.dual();
1114      CoinFillN(rowObjective,numberRows,0.0);
1115      CoinDisjointCopyN(m.getObjCoefficients(),numberColumns,objective);
1116      solution.dual();
1117    }
1118  }
1119  // test unbounded
1120  {   
1121    CoinMpsIO m;
1122    std::string fn = mpsDir+"brandy";
1123    m.readMps(fn.c_str(),"mps");
1124    ClpSimplex solution;
1125    // do twice - without and with scaling
1126    int iPass;
1127    for (iPass=0;iPass<2;iPass++) {
1128      solution.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
1129                      m.getObjCoefficients(),
1130                      m.getRowLower(),m.getRowUpper());
1131      if (iPass)
1132        solution.scaling();
1133      solution.setOptimizationDirection(-1);
1134      // test unbounded and ray
1135#ifdef DUAL
1136      solution.setDualBound(100.0);
1137      solution.dual();
1138#else
1139      solution.primal();
1140#endif
1141      assert (solution.status()==2);
1142      int numberColumns = solution.numberColumns();
1143      int numberRows = solution.numberRows();
1144      double * lower = solution.columnLower();
1145      double * upper = solution.columnUpper();
1146      double * sol = solution.primalColumnSolution();
1147      double * ray = solution.unboundedRay();
1148      double * objective = solution.objective();
1149      double objChange=0.0;
1150      int iRow , iColumn;
1151      // make sure feasible and columns form ray
1152      for (iColumn=0;iColumn<numberColumns;iColumn++) {
1153        double value = sol[iColumn];
1154        assert(value<upper[iColumn]+1.0e-8);
1155        assert(value>lower[iColumn]-1.0e-8);
1156        value = ray[iColumn];
1157        if (value>0.0)
1158          assert(upper[iColumn]>1.0e30);
1159        else if (value<0.0)
1160          assert(lower[iColumn]<-1.0e30);
1161        objChange += value*objective[iColumn];
1162      }
1163      // make sure increasing objective
1164      assert(objChange>0.0);
1165      double * result = new double[numberRows];
1166      CoinFillN ( result, numberRows,0.0);
1167      solution.matrix()->times(sol, result);
1168      lower = solution.rowLower();
1169      upper = solution.rowUpper();
1170      sol = solution.primalRowSolution();
1171      for (iRow=0;iRow<numberRows;iRow++) {
1172        double value = result[iRow];
1173        assert(eq(value,sol[iRow]));
1174        assert(value<upper[iRow]+2.0e-8);
1175        assert(value>lower[iRow]-2.0e-8);
1176      }
1177      CoinFillN ( result, numberRows,0.0);
1178      solution.matrix()->times(ray, result);
1179      // there may be small differences (especially if scaled)
1180      for (iRow=0;iRow<numberRows;iRow++) {
1181        double value = result[iRow];
1182        if (value>1.0e-8)
1183          assert(upper[iRow]>1.0e30);
1184        else if (value<-1.0e-8)
1185          assert(lower[iRow]<-1.0e30);
1186      }
1187      delete [] result;
1188      delete [] ray;
1189    }
1190  }
1191  // test infeasible
1192  {   
1193    CoinMpsIO m;
1194    std::string fn = mpsDir+"brandy";
1195    m.readMps(fn.c_str(),"mps");
1196    ClpSimplex solution;
1197    // do twice - without and with scaling
1198    int iPass;
1199    for (iPass=0;iPass<2;iPass++) {
1200      solution.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
1201                      m.getObjCoefficients(),
1202                      m.getRowLower(),m.getRowUpper());
1203      if (iPass)
1204        solution.scaling();
1205      // test infeasible and ray
1206      solution.columnUpper()[0]=0.0;
1207#ifdef DUAL
1208      solution.setDualBound(100.0);
1209      solution.dual();
1210#else
1211      solution.primal();
1212#endif
1213      assert (solution.status()==1);
1214      int numberColumns = solution.numberColumns();
1215      int numberRows = solution.numberRows();
1216      double * lower = solution.rowLower();
1217      double * upper = solution.rowUpper();
1218      double * ray = solution.infeasibilityRay();
1219      assert(ray);
1220      // construct proof of infeasibility
1221      int iRow , iColumn;
1222      double lo=0.0,up=0.0;
1223      int nl=0,nu=0;
1224      for (iRow=0;iRow<numberRows;iRow++) {
1225        if (lower[iRow]>-1.0e20) {
1226          lo += ray[iRow]*lower[iRow];
1227        } else {
1228          if (ray[iRow]>1.0e-8) 
1229            nl++;
1230        }
1231        if (upper[iRow]<1.0e20) {
1232          up += ray[iRow]*upper[iRow];
1233        } else {
1234          if (ray[iRow]>1.0e-8) 
1235            nu++;
1236        }
1237      }
1238      if (nl)
1239        lo=-1.0e100;
1240      if (nu)
1241        up=1.0e100;
1242      double * result = new double[numberColumns];
1243      double lo2=0.0,up2=0.0;
1244      CoinFillN ( result, numberColumns,0.0);
1245      solution.matrix()->transposeTimes(ray, result);
1246      lower = solution.columnLower();
1247      upper = solution.columnUpper();
1248      nl=nu=0;
1249      for (iColumn=0;iColumn<numberColumns;iColumn++) {
1250        if (result[iColumn]>1.0e-8) {
1251          if (lower[iColumn]>-1.0e20)
1252            lo2 += result[iColumn]*lower[iColumn];
1253          else
1254            nl++;
1255          if (upper[iColumn]<1.0e20)
1256            up2 += result[iColumn]*upper[iColumn];
1257          else
1258            nu++;
1259        } else if (result[iColumn]<-1.0e-8) {
1260          if (lower[iColumn]>-1.0e20)
1261            up2 += result[iColumn]*lower[iColumn];
1262          else
1263            nu++;
1264          if (upper[iColumn]<1.0e20)
1265            lo2 += result[iColumn]*upper[iColumn];
1266          else
1267            nl++;
1268        }
1269      }
1270      if (nl)
1271        lo2=-1.0e100;
1272      if (nu)
1273        up2=1.0e100;
1274      // make sure inconsistency
1275      assert(lo2>up||up2<lo);
1276      delete [] result;
1277      delete [] ray;
1278    }
1279  }
1280  // test delete and add
1281  {   
1282    CoinMpsIO m;
1283    std::string fn = mpsDir+"brandy";
1284    m.readMps(fn.c_str(),"mps");
1285    ClpSimplex solution;
1286    solution.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
1287                         m.getObjCoefficients(),
1288                         m.getRowLower(),m.getRowUpper());
1289    solution.dual();
1290    CoinRelFltEq eq(1.0e-8);
1291    assert(eq(solution.objectiveValue(),1.5185098965e+03));
1292
1293    int numberColumns = solution.numberColumns();
1294    int numberRows = solution.numberRows();
1295    double * saveObj = new double [numberColumns];
1296    double * saveLower = new double[numberRows+numberColumns];
1297    double * saveUpper = new double[numberRows+numberColumns];
1298    int * which = new int [numberRows+numberColumns];
1299
1300    int numberElements = m.getMatrixByCol()->getNumElements();
1301    int * starts = new int[numberRows+numberColumns];
1302    int * index = new int[numberElements];
1303    double * element = new double[numberElements];
1304
1305    const CoinBigIndex * startM;
1306    const int * lengthM;
1307    const int * indexM;
1308    const double * elementM;
1309
1310    int n,nel;
1311
1312    // delete non basic columns
1313    n=0;
1314    nel=0;
1315    int iRow , iColumn;
1316    const double * lower = m.getColLower();
1317    const double * upper = m.getColUpper();
1318    const double * objective = m.getObjCoefficients();
1319    startM = m.getMatrixByCol()->getVectorStarts();
1320    lengthM = m.getMatrixByCol()->getVectorLengths();
1321    indexM = m.getMatrixByCol()->getIndices();
1322    elementM = m.getMatrixByCol()->getElements();
1323    starts[0]=0;
1324    for (iColumn=0;iColumn<numberColumns;iColumn++) {
1325      if (solution.getColumnStatus(iColumn)!=ClpSimplex::basic) {
1326        saveObj[n]=objective[iColumn];
1327        saveLower[n]=lower[iColumn];
1328        saveUpper[n]=upper[iColumn];
1329        int j;
1330        for (j=startM[iColumn];j<startM[iColumn]+lengthM[iColumn];j++) {
1331          index[nel]=indexM[j];
1332          element[nel++]=elementM[j];
1333        }
1334        which[n++]=iColumn;
1335        starts[n]=nel;
1336      }
1337    }
1338    solution.deleteColumns(n,which);
1339    solution.dual();
1340    // Put back
1341    solution.addColumns(n,saveLower,saveUpper,saveObj,
1342                        starts,index,element);
1343    solution.dual();
1344    assert(eq(solution.objectiveValue(),1.5185098965e+03));
1345    // Delete all columns and add back
1346    n=0;
1347    nel=0;
1348    starts[0]=0;
1349    lower = m.getColLower();
1350    upper = m.getColUpper();
1351    objective = m.getObjCoefficients();
1352    for (iColumn=0;iColumn<numberColumns;iColumn++) {
1353      saveObj[n]=objective[iColumn];
1354      saveLower[n]=lower[iColumn];
1355      saveUpper[n]=upper[iColumn];
1356      int j;
1357      for (j=startM[iColumn];j<startM[iColumn]+lengthM[iColumn];j++) {
1358        index[nel]=indexM[j];
1359        element[nel++]=elementM[j];
1360      }
1361      which[n++]=iColumn;
1362      starts[n]=nel;
1363    }
1364    solution.deleteColumns(n,which);
1365    solution.dual();
1366    // Put back
1367    solution.addColumns(n,saveLower,saveUpper,saveObj,
1368                        starts,index,element);
1369    solution.dual();
1370    assert(eq(solution.objectiveValue(),1.5185098965e+03));
1371
1372    // reload with original
1373    solution.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
1374                         m.getObjCoefficients(),
1375                         m.getRowLower(),m.getRowUpper());
1376    // delete half rows
1377    n=0;
1378    nel=0;
1379    lower = m.getRowLower();
1380    upper = m.getRowUpper();
1381    startM = m.getMatrixByRow()->getVectorStarts();
1382    lengthM = m.getMatrixByRow()->getVectorLengths();
1383    indexM = m.getMatrixByRow()->getIndices();
1384    elementM = m.getMatrixByRow()->getElements();
1385    starts[0]=0;
1386    for (iRow=0;iRow<numberRows;iRow++) {
1387      if ((iRow&1)==0) {
1388        saveLower[n]=lower[iRow];
1389        saveUpper[n]=upper[iRow];
1390        int j;
1391        for (j=startM[iRow];j<startM[iRow]+lengthM[iRow];j++) {
1392          index[nel]=indexM[j];
1393          element[nel++]=elementM[j];
1394        }
1395        which[n++]=iRow;
1396        starts[n]=nel;
1397      }
1398    }
1399    solution.deleteRows(n,which);
1400    solution.dual();
1401    // Put back
1402    solution.addRows(n,saveLower,saveUpper,
1403                        starts,index,element);
1404    solution.dual();
1405    assert(eq(solution.objectiveValue(),1.5185098965e+03));
1406    solution.writeMps("yy.mps");
1407    // Delete all rows
1408    n=0;
1409    nel=0;
1410    lower = m.getRowLower();
1411    upper = m.getRowUpper();
1412    starts[0]=0;
1413    for (iRow=0;iRow<numberRows;iRow++) {
1414      saveLower[n]=lower[iRow];
1415      saveUpper[n]=upper[iRow];
1416      int j;
1417      for (j=startM[iRow];j<startM[iRow]+lengthM[iRow];j++) {
1418        index[nel]=indexM[j];
1419        element[nel++]=elementM[j];
1420      }
1421      which[n++]=iRow;
1422      starts[n]=nel;
1423    }
1424    solution.deleteRows(n,which);
1425    solution.dual();
1426    // Put back
1427    solution.addRows(n,saveLower,saveUpper,
1428                        starts,index,element);
1429    solution.dual();
1430    solution.writeMps("xx.mps");
1431    assert(eq(solution.objectiveValue(),1.5185098965e+03));
1432    // Zero out status array to give some interest
1433    memset(solution.statusArray()+numberColumns,0,numberRows);
1434    solution.primal(1);
1435    assert(eq(solution.objectiveValue(),1.5185098965e+03));
1436    // Delete all columns and rows
1437    n=0;
1438    for (iColumn=0;iColumn<numberColumns;iColumn++) {
1439      which[n++]=iColumn;
1440      starts[n]=nel;
1441    }
1442    solution.deleteColumns(n,which);
1443    n=0;
1444    for (iRow=0;iRow<numberRows;iRow++) {
1445      which[n++]=iRow;
1446      starts[n]=nel;
1447    }
1448    solution.deleteRows(n,which);
1449
1450    delete [] saveObj;
1451    delete [] saveLower;
1452    delete [] saveUpper;
1453    delete [] which;
1454    delete [] starts;
1455    delete [] index;
1456    delete [] element;
1457  }
1458#if 1
1459  // Test barrier
1460  {
1461    CoinMpsIO m;
1462    std::string fn = mpsDir+"exmip1";
1463    m.readMps(fn.c_str(),"mps");
1464    ClpInterior solution;
1465    solution.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
1466                         m.getObjCoefficients(),
1467                         m.getRowLower(),m.getRowUpper());
1468    solution.primalDual();
1469  }
1470#endif
1471  // test network
1472#define QUADRATIC
1473  if (1) {   
1474    std::string fn = mpsDir+"input.130";
1475    int numberColumns;
1476    int numberRows;
1477   
1478    FILE * fp = fopen(fn.c_str(),"r");
1479    if (!fp) {
1480      // Try in Data/Sample
1481      fn = "Data/Sample/input.130";
1482      fp = fopen(fn.c_str(),"r");
1483    }
1484    if (!fp) {
1485      fprintf(stderr,"Unable to open file input.130 in mpsDir or Data/Sample directory\n");
1486    } else {
1487      int problem;
1488      char temp[100];
1489      // read and skip
1490      fscanf(fp,"%s",temp);
1491      assert (!strcmp(temp,"BEGIN"));
1492      fscanf(fp,"%*s %*s %d %d %*s %*s %d %*s",&problem, &numberRows, 
1493             &numberColumns);
1494      // scan down to SUPPLY
1495      while (fgets(temp,100,fp)) {
1496        if (!strncmp(temp,"SUPPLY",6))
1497          break;
1498      }
1499      if (strncmp(temp,"SUPPLY",6)) {
1500        fprintf(stderr,"Unable to find SUPPLY\n");
1501        exit(2);
1502      }
1503      // get space for rhs
1504      double * lower = new double[numberRows];
1505      double * upper = new double[numberRows];
1506      int i;
1507      for (i=0;i<numberRows;i++) {
1508        lower[i]=0.0;
1509        upper[i]=0.0;
1510      }
1511      // ***** Remember to convert to C notation
1512      while (fgets(temp,100,fp)) {
1513        int row;
1514        int value;
1515        if (!strncmp(temp,"ARCS",4))
1516          break;
1517        sscanf(temp,"%d %d",&row,&value);
1518        upper[row-1]=-value;
1519        lower[row-1]=-value;
1520      }
1521      if (strncmp(temp,"ARCS",4)) {
1522        fprintf(stderr,"Unable to find ARCS\n");
1523        exit(2);
1524      }
1525      // number of columns may be underestimate
1526      int * head = new int[2*numberColumns];
1527      int * tail = new int[2*numberColumns];
1528      int * ub = new int[2*numberColumns];
1529      int * cost = new int[2*numberColumns];
1530      // ***** Remember to convert to C notation
1531      numberColumns=0;
1532      while (fgets(temp,100,fp)) {
1533        int iHead;
1534        int iTail;
1535        int iUb;
1536        int iCost;
1537        if (!strncmp(temp,"DEMAND",6))
1538          break;
1539        sscanf(temp,"%d %d %d %d",&iHead,&iTail,&iCost,&iUb);
1540        iHead--;
1541        iTail--;
1542        head[numberColumns]=iHead;
1543        tail[numberColumns]=iTail;
1544        ub[numberColumns]=iUb;
1545        cost[numberColumns]=iCost;
1546        numberColumns++;
1547      }
1548      if (strncmp(temp,"DEMAND",6)) {
1549        fprintf(stderr,"Unable to find DEMAND\n");
1550        exit(2);
1551      }
1552      // ***** Remember to convert to C notation
1553      while (fgets(temp,100,fp)) {
1554        int row;
1555        int value;
1556        if (!strncmp(temp,"END",3))
1557          break;
1558        sscanf(temp,"%d %d",&row,&value);
1559        upper[row-1]=value;
1560        lower[row-1]=value;
1561      }
1562      if (strncmp(temp,"END",3)) {
1563        fprintf(stderr,"Unable to find END\n");
1564        exit(2);
1565      }
1566      printf("Problem %d has %d rows and %d columns\n",problem,
1567             numberRows,numberColumns);
1568      fclose(fp);
1569      ClpSimplex  model;
1570      // now build model
1571     
1572      double * objective =new double[numberColumns];
1573      double * lowerColumn = new double[numberColumns];
1574      double * upperColumn = new double[numberColumns];
1575     
1576      double * element = new double [2*numberColumns];
1577      CoinBigIndex * start = new CoinBigIndex [numberColumns+1];
1578      int * row = new int[2*numberColumns];
1579      start[numberColumns]=2*numberColumns;
1580      for (i=0;i<numberColumns;i++) {
1581        start[i]=2*i;
1582        element[2*i]=-1.0;
1583        element[2*i+1]=1.0;
1584        row[2*i]=head[i];
1585        row[2*i+1]=tail[i];
1586        lowerColumn[i]=0.0;
1587        upperColumn[i]=ub[i];
1588        objective[i]=cost[i];
1589      }
1590      // Create Packed Matrix
1591      CoinPackedMatrix matrix;
1592      int * lengths = NULL;
1593      matrix.assignMatrix(true,numberRows,numberColumns,
1594                          2*numberColumns,element,row,start,lengths);
1595      // load model
1596      model.loadProblem(matrix,
1597                        lowerColumn,upperColumn,objective,
1598                        lower,upper);
1599      model.factorization()->maximumPivots(200+model.numberRows()/100);
1600      model.createStatus();
1601      double time1 = CoinCpuTime();
1602      model.dual();
1603      std::cout<<"Network problem, ClpPackedMatrix took "<<CoinCpuTime()-time1<<" seconds"<<std::endl;
1604      ClpPlusMinusOneMatrix * plusMinus = new ClpPlusMinusOneMatrix(matrix);
1605      assert (plusMinus->getIndices()); // would be zero if not +- one
1606      //ClpPlusMinusOneMatrix *plusminus_matrix;
1607
1608      //plusminus_matrix = new ClpPlusMinusOneMatrix;
1609
1610      //plusminus_matrix->passInCopy(numberRows, numberColumns, true, plusMinus->getMutableIndices(),
1611      //                         plusMinus->startPositive(),plusMinus->startNegative());
1612      model.loadProblem(*plusMinus,
1613                        lowerColumn,upperColumn,objective,
1614                        lower,upper);
1615      //model.replaceMatrix( plusminus_matrix , true);
1616      delete plusMinus;
1617      //model.createStatus();
1618      //model.initialSolve();
1619      //model.writeMps("xx.mps");
1620     
1621      model.factorization()->maximumPivots(200+model.numberRows()/100);
1622      model.createStatus();
1623      time1 = CoinCpuTime();
1624      model.dual();
1625      std::cout<<"Network problem, ClpPlusMinusOneMatrix took "<<CoinCpuTime()-time1<<" seconds"<<std::endl;
1626      ClpNetworkMatrix network(numberColumns,head,tail);
1627      model.loadProblem(network,
1628                        lowerColumn,upperColumn,objective,
1629                        lower,upper);
1630     
1631      model.factorization()->maximumPivots(200+model.numberRows()/100);
1632      model.createStatus();
1633      time1 = CoinCpuTime();
1634      model.dual();
1635      std::cout<<"Network problem, ClpNetworkMatrix took "<<CoinCpuTime()-time1<<" seconds"<<std::endl;
1636      delete [] lower;
1637      delete [] upper;
1638      delete [] head;
1639      delete [] tail;
1640      delete [] ub;
1641      delete [] cost;
1642      delete [] objective;
1643      delete [] lowerColumn;
1644      delete [] upperColumn;
1645    }
1646  }
1647#ifdef QUADRATIC
1648  // Test quadratic to solve linear
1649  if (1) {   
1650    CoinMpsIO m;
1651    std::string fn = mpsDir+"exmip1";
1652    m.readMps(fn.c_str(),"mps");
1653    ClpSimplex solution;
1654    solution.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
1655                         m.getObjCoefficients(),
1656                         m.getRowLower(),m.getRowUpper());
1657    //solution.dual();
1658    // get quadratic part
1659    int numberColumns=solution.numberColumns();
1660    int * start=new int [numberColumns+1];
1661    int * column = new int[numberColumns];
1662    double * element = new double[numberColumns];
1663    int i;
1664    start[0]=0;
1665    int n=0;
1666    int kk=numberColumns-1;
1667    int kk2=numberColumns-1;
1668    for (i=0;i<numberColumns;i++) {
1669      if (i>=kk) {
1670        column[n]=i;
1671        if (i>=kk2)
1672          element[n]=1.0e-1;
1673        else
1674          element[n]=0.0;
1675        n++;
1676      }
1677      start[i+1]=n;
1678    }
1679    // Load up objective
1680    solution.loadQuadraticObjective(numberColumns,start,column,element);
1681    delete [] start;
1682    delete [] column;
1683    delete [] element;
1684    //solution.quadraticSLP(50,1.0e-4);
1685    CoinRelFltEq eq(1.0e-4);
1686    //assert(eq(objValue,820.0));
1687    //solution.setLogLevel(63);
1688    solution.primal();
1689    printSol(solution);
1690    //assert(eq(objValue,3.2368421));
1691    //exit(77);
1692  }
1693  // Test quadratic
1694  if (1) {   
1695    CoinMpsIO m;
1696    std::string fn = mpsDir+"share2qp";
1697    //fn = "share2qpb";
1698    m.readMps(fn.c_str(),"mps");
1699    ClpSimplex model;
1700    model.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
1701                         m.getObjCoefficients(),
1702                         m.getRowLower(),m.getRowUpper());
1703    model.dual();
1704    // get quadratic part
1705    int * start=NULL;
1706    int * column = NULL;
1707    double * element = NULL;
1708    m.readQuadraticMps(NULL,start,column,element,2);
1709    int column2[200];
1710    double element2[200];
1711    int start2[80];
1712    int j;
1713    start2[0]=0;
1714    int nel=0;
1715    bool good=false;
1716    for (j=0;j<79;j++) {
1717      if (start[j]==start[j+1]) {
1718        column2[nel]=j;
1719        element2[nel]=0.0;
1720        nel++;
1721      } else {
1722        int i;
1723        for (i=start[j];i<start[j+1];i++) {
1724          column2[nel]=column[i];
1725          element2[nel++]=element[i];
1726        }
1727      }
1728      start2[j+1]=nel;
1729    }
1730    // Load up objective
1731    if (good)
1732      model.loadQuadraticObjective(model.numberColumns(),start2,column2,element2);
1733    else
1734      model.loadQuadraticObjective(model.numberColumns(),start,column,element);
1735    delete [] start;
1736    delete [] column;
1737    delete [] element;
1738    int numberColumns=model.numberColumns();
1739    model.scaling(0);
1740#if 0
1741    model.nonlinearSLP(50,1.0e-4);
1742#else
1743    // Get feasible
1744    ClpObjective * saveObjective = model.objectiveAsObject()->clone();
1745    ClpLinearObjective zeroObjective(NULL,numberColumns);
1746    model.setObjective(&zeroObjective);
1747    model.dual();
1748    model.setObjective(saveObjective);
1749    delete saveObjective;
1750#endif
1751    //model.setLogLevel(63);
1752    //exit(77);
1753    model.setFactorizationFrequency(10);
1754    model.primal();
1755    printSol(model);
1756    double objValue = model.getObjValue();
1757    CoinRelFltEq eq(1.0e-4);
1758    assert(eq(objValue,-400.92));
1759    // and again for barrier
1760    model.barrier(false);
1761    //printSol(model);
1762    model.allSlackBasis();
1763    model.primal();
1764    //printSol(model);
1765  }
1766  if (0) {   
1767    CoinMpsIO m;
1768    std::string fn = "./beale";
1769    //fn = "./jensen";
1770    m.readMps(fn.c_str(),"mps");
1771    ClpSimplex solution;
1772    solution.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
1773                         m.getObjCoefficients(),
1774                         m.getRowLower(),m.getRowUpper());
1775    solution.setDblParam(ClpObjOffset,m.objectiveOffset());
1776    solution.dual();
1777    // get quadratic part
1778    int * start=NULL;
1779    int * column = NULL;
1780    double * element = NULL;
1781    m.readQuadraticMps(NULL,start,column,element,2);
1782    // Load up objective
1783    solution.loadQuadraticObjective(solution.numberColumns(),start,column,element);
1784    delete [] start;
1785    delete [] column;
1786    delete [] element;
1787    solution.primal(1);
1788    solution.nonlinearSLP(50,1.0e-4);
1789    double objValue = solution.getObjValue();
1790    CoinRelFltEq eq(1.0e-4);
1791    assert(eq(objValue,0.5));
1792    solution.primal();
1793    objValue = solution.getObjValue();
1794    assert(eq(objValue,0.5));
1795  }
1796#endif 
1797}
Note: See TracBrowser for help on using the repository browser.