source: stable/1.6/Clp/src/unitTest.cpp @ 1609

Last change on this file since 1609 was 1103, checked in by forrest, 13 years ago

scale binvrow

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