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

Last change on this file since 1071 was 1056, checked in by forrest, 13 years ago

out ambiguities

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