source: trunk/Clp/src/unitTest.cpp @ 1321

Last change on this file since 1321 was 1321, checked in by forrest, 11 years ago

out compiler warnings and stability improvements

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