source: trunk/Clp/test/unitTest.cpp @ 754

Last change on this file since 754 was 754, checked in by andreasw, 14 years ago

first version

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