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

Last change on this file since 1361 was 1361, checked in by forrest, 12 years ago

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