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

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

minor changes

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 75.3 KB
Line 
1/* $Id: unitTest.cpp 1434 2009-09-24 20:48:01Z forrest $ */
2// Copyright (C) 2002, International Business Machines
3// Corporation and others.  All Rights Reserved.
4
5#ifdef NDEBUG
6#undef NDEBUG
7#endif
8
9#include "ClpConfig.h"
10#include "CoinPragma.hpp"
11#include <cassert>
12#include <cstdio>
13#include <cstdlib>
14#include <cmath>
15#include <cfloat>
16#include <string>
17#include <iostream>
18
19#include "CoinMpsIO.hpp"
20#include "CoinPackedMatrix.hpp"
21#include "CoinPackedVector.hpp"
22#include "CoinStructuredModel.hpp"
23#include "CoinHelperFunctions.hpp"
24#include "CoinTime.hpp"
25
26#include "ClpFactorization.hpp"
27#include "ClpSimplex.hpp"
28#include "ClpSimplexOther.hpp"
29#include "ClpSimplexNonlinear.hpp"
30#include "ClpInterior.hpp"
31#include "ClpLinearObjective.hpp"
32#include "ClpDualRowSteepest.hpp"
33#include "ClpDualRowDantzig.hpp"
34#include "ClpPrimalColumnSteepest.hpp"
35#include "ClpPrimalColumnDantzig.hpp"
36#include "ClpParameters.hpp"
37#include "ClpNetworkMatrix.hpp"
38#include "ClpPlusMinusOneMatrix.hpp"
39#include "MyMessageHandler.hpp"
40#include "MyEventHandler.hpp"
41
42#include "ClpPresolve.hpp"
43#include "Idiot.hpp"
44
45
46//#############################################################################
47
48
49// Function Prototypes. Function definitions is in this file.
50void testingMessage( const char * const msg );
51#if UFL_BARRIER
52static int barrierAvailable=1;
53static std::string nameBarrier="barrier-UFL";
54#elif WSSMP_BARRIER
55static int barrierAvailable=2;
56static std::string nameBarrier="barrier-WSSMP";
57#elif MUMPS_BARRIER
58static int barrierAvailable=3;
59static std::string nameBarrier="barrier-MUMPS";
60#else
61static int barrierAvailable=0;
62static std::string nameBarrier="barrier-slow";
63#endif
64#define NUMBER_ALGORITHMS 12
65// If you just want a subset then set some to 1
66static int switchOff[NUMBER_ALGORITHMS]={0,0,0,0,0,0,0,0,0,0,0,0};
67// shortName - 0 no , 1 yes
68ClpSolve setupForSolve(int algorithm, std::string & nameAlgorithm,
69                       int shortName)
70{
71  ClpSolve solveOptions;
72  /* algorithms are
73     0 barrier
74     1 dual with volumne crash
75     2,3 dual with and without crash
76     4,5 primal with and without
77     6,7 automatic with and without
78     8,9 primal with idiot 1 and 5
79     10,11 primal with 70, dual with volume
80  */
81  switch (algorithm) {
82  case 0:
83    if (shortName)
84      nameAlgorithm="ba";
85    else
86      nameAlgorithm="nameBarrier";
87    solveOptions.setSolveType(ClpSolve::useBarrier);
88    if (barrierAvailable==1) 
89      solveOptions.setSpecialOption(4,4);
90    else if (barrierAvailable==2) 
91      solveOptions.setSpecialOption(4,2);
92    break;
93  case 1:
94#ifdef COIN_HAS_VOL
95    if (shortName)
96      nameAlgorithm="du-vol-50";
97    else
98      nameAlgorithm="dual-volume-50";
99    solveOptions.setSolveType(ClpSolve::useDual);
100    solveOptions.setSpecialOption(0,2,50); // volume
101#else 
102      solveOptions.setSolveType(ClpSolve::notImplemented);
103#endif
104    break;
105  case 2:
106    if (shortName)
107      nameAlgorithm="du-cr";
108    else
109      nameAlgorithm="dual-crash";
110    solveOptions.setSolveType(ClpSolve::useDual);
111    solveOptions.setSpecialOption(0,1);
112    break;
113  case 3:
114    if (shortName)
115      nameAlgorithm="du";
116    else
117      nameAlgorithm="dual";
118    solveOptions.setSolveType(ClpSolve::useDual);
119    break;
120  case 4:
121    if (shortName)
122      nameAlgorithm="pr-cr";
123    else
124      nameAlgorithm="primal-crash";
125    solveOptions.setSolveType(ClpSolve::usePrimal);
126    solveOptions.setSpecialOption(1,1);
127    break;
128  case 5:
129    if (shortName)
130      nameAlgorithm="pr";
131    else
132      nameAlgorithm="primal";
133    solveOptions.setSolveType(ClpSolve::usePrimal);
134    break;
135  case 6:
136    if (shortName)
137      nameAlgorithm="au-cr";
138    else
139      nameAlgorithm="either-crash";
140    solveOptions.setSolveType(ClpSolve::automatic);
141    solveOptions.setSpecialOption(1,1);
142    break;
143  case 7:
144    if (shortName)
145      nameAlgorithm="au";
146    else
147      nameAlgorithm="either";
148    solveOptions.setSolveType(ClpSolve::automatic);
149    break;
150  case 8:
151    if (shortName)
152      nameAlgorithm="pr-id-1";
153    else
154      nameAlgorithm="primal-idiot-1";
155    solveOptions.setSolveType(ClpSolve::usePrimalorSprint);
156    solveOptions.setSpecialOption(1,2,1); // idiot
157    break;
158  case 9:
159    if (shortName)
160      nameAlgorithm="pr-id-5";
161    else
162      nameAlgorithm="primal-idiot-5";
163    solveOptions.setSolveType(ClpSolve::usePrimalorSprint);
164    solveOptions.setSpecialOption(1,2,5); // idiot
165    break;
166  case 10:
167    if (shortName)
168      nameAlgorithm="pr-id-70";
169    else
170      nameAlgorithm="primal-idiot-70";
171    solveOptions.setSolveType(ClpSolve::usePrimalorSprint);
172    solveOptions.setSpecialOption(1,2,70); // idiot
173    break;
174  case 11:
175#ifdef COIN_HAS_VOL
176    if (shortName)
177      nameAlgorithm="du-vol";
178    else
179      nameAlgorithm="dual-volume";
180    solveOptions.setSolveType(ClpSolve::useDual);
181    solveOptions.setSpecialOption(0,2,3000); // volume
182#else 
183      solveOptions.setSolveType(ClpSolve::notImplemented);
184#endif
185    break;
186  default:
187    abort();
188  }
189  if (shortName) {
190    // can switch off
191    if (switchOff[algorithm])
192      solveOptions.setSolveType(ClpSolve::notImplemented);
193  }
194  return solveOptions;
195}
196static void printSol(ClpSimplex & model) 
197{
198  int numberRows = model.numberRows();
199  int numberColumns = model.numberColumns();
200 
201  double * rowPrimal = model.primalRowSolution();
202  double * rowDual = model.dualRowSolution();
203  double * rowLower = model.rowLower();
204  double * rowUpper = model.rowUpper();
205 
206  int iRow;
207  double objValue = model.getObjValue();
208  printf("Objvalue %g Rows (%d)\n",objValue,numberRows);
209  for (iRow=0;iRow<numberRows;iRow++) {
210    printf("%d primal %g dual %g low %g up %g\n",
211           iRow,rowPrimal[iRow],rowDual[iRow],
212           rowLower[iRow],rowUpper[iRow]);
213  }
214  double * columnPrimal = model.primalColumnSolution();
215  double * columnDual = model.dualColumnSolution();
216  double * columnLower = model.columnLower();
217  double * columnUpper = model.columnUpper();
218  double offset;
219  //const double * gradient = model.objectiveAsObject()->gradient(&model,
220  //                                                       columnPrimal,offset,true,1);
221  const double * gradient = model.objective(columnPrimal,offset);
222  int iColumn;
223  objValue = -offset-model.objectiveOffset();
224  printf("offset %g (%g)\n",offset,model.objectiveOffset());
225  printf("Columns (%d)\n",numberColumns);
226  for (iColumn=0;iColumn<numberColumns;iColumn++) {
227    printf("%d primal %g dual %g low %g up %g\n",
228           iColumn,columnPrimal[iColumn],columnDual[iColumn],
229           columnLower[iColumn],columnUpper[iColumn]);
230    objValue += columnPrimal[iColumn]*gradient[iColumn];
231    if (fabs(columnPrimal[iColumn]*gradient[iColumn])>1.0e-8)
232      printf("obj -> %g gradient %g\n",objValue,gradient[iColumn]);
233  }
234  printf("Computed objective %g\n",objValue);
235}
236
237void usage(const std::string& key)
238{
239  std::cerr
240    <<"Undefined parameter \"" <<key <<"\".\n"
241    <<"Correct usage: \n"
242    <<"  clp [-dirSample=V1] [-dirNetlib=V2] [-netlib]\n"
243    <<"  where:\n"
244    <<"    -dirSample: directory containing mps test files\n"
245    <<"        Default value V1=\"../../Data/Sample\"\n"
246    <<"    -dirNetlib: directory containing netlib files\"\n"
247    <<"        Default value V2=\"../../Data/Netlib\"\n"
248    <<"    -netlib\n"
249    <<"        If specified, then netlib testset run as well as the nitTest.\n";
250}
251
252//----------------------------------------------------------------
253int mainTest (int argc, const char *argv[],int algorithm,
254              ClpSimplex empty, ClpSolve solveOptionsIn,
255              int switchOffValue,bool doVector)
256{
257  int i;
258
259  if (switchOffValue>0) {
260    // switch off some
261    int iTest;
262    for (iTest=0;iTest<NUMBER_ALGORITHMS;iTest++) {
263#ifndef NDEBUG
264      int bottom = switchOffValue%10;
265      assert (bottom==0||bottom==1);
266#endif
267      switchOffValue/= 10;
268      switchOff[iTest]=0;
269    }
270  }
271
272  // define valid parameter keywords
273  std::set<std::string> definedKeyWords;
274  definedKeyWords.insert("-dirSample");
275  definedKeyWords.insert("-dirNetlib");
276  definedKeyWords.insert("-netlib");
277
278  // Create a map of parameter keys and associated data
279  std::map<std::string,std::string> parms;
280  for ( i=1; i<argc; i++ ) {
281    std::string parm(argv[i]);
282    std::string key,value;
283    std::string::size_type  eqPos = parm.find('=');
284
285    // Does parm contain and '='
286    if ( eqPos==std::string::npos ) {
287      //Parm does not contain '='
288      key = parm;
289    }
290    else {
291      key=parm.substr(0,eqPos);
292      value=parm.substr(eqPos+1);
293    }
294
295    // Is specifed key valid?
296    if ( definedKeyWords.find(key) == definedKeyWords.end() ) {
297      // invalid key word.
298      // Write help text
299      usage(key);
300      return 1;
301    }
302    parms[key]=value;
303  }
304 
305  const char dirsep =  CoinFindDirSeparator();
306  // Set directory containing mps data files.
307  std::string dirSample;
308  if (parms.find("-dirSample") != parms.end())
309    dirSample=parms["-dirSample"];
310  else 
311    dirSample = dirsep == '/' ? "../../Data/Sample/" : "..\\..\\Data\\Sample\\";
312 
313  // Set directory containing netlib data files.
314  std::string dirNetlib;
315  if (parms.find("-dirNetlib") != parms.end())
316    dirNetlib=parms["-dirNetlib"];
317  else 
318    dirNetlib = dirsep == '/' ? "../../Data/Netlib/" : "..\\..\\Data\\Netlib\\";
319  if (!empty.numberRows()) {
320    testingMessage( "Testing ClpSimplex\n" );
321    ClpSimplexUnitTest(dirSample);
322  }
323  if (parms.find("-netlib") != parms.end()||empty.numberRows())
324  {
325    unsigned int m;
326   
327    // Define test problems:
328    //   mps names,
329    //   maximization or minimization,
330    //   Number of rows and columns in problem, and
331    //   objective function value
332    std::vector<std::string> mpsName;
333    std::vector<bool> min;
334    std::vector<int> nRows;
335    std::vector<int> nCols;
336    std::vector<double> objValue;
337    std::vector<double> objValueTol;
338    // 100 added means no presolve
339    std::vector<int> bestStrategy;
340    if(empty.numberRows()) {
341      std::string alg;
342      for (int iTest=0;iTest<NUMBER_ALGORITHMS;iTest++) {
343        ClpSolve solveOptions=setupForSolve(iTest,alg,0);
344        printf("%d %s ",iTest,alg.c_str());
345        if (switchOff[iTest]) 
346          printf("skipped by user\n");
347        else if(solveOptions.getSolveType()==ClpSolve::notImplemented)
348          printf("skipped as not available\n");
349        else
350          printf("will be tested\n");
351      }
352    }
353    if (!empty.numberRows()) {
354      mpsName.push_back("25fv47");
355      min.push_back(true);
356      nRows.push_back(822);
357      nCols.push_back(1571);
358      objValueTol.push_back(1.E-10);
359      objValue.push_back(5.5018458883E+03);
360      bestStrategy.push_back(0);
361      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);
362      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);
363      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);
364      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);
365      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);
366      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);
367      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);
368      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);
369      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);
370      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);
371      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);
372      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);
373      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);
374      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);
375      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);
376      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);
377      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);
378      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);
379      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);
380      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);
381      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);
382      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);
383      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);
384      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);
385      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);
386      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);
387      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);
388      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);
389      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.
390      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);
391      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);
392      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);
393      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);
394      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);
395      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);
396      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);
397      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);
398      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);
399      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);
400      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);
401      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);
402      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);
403      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);
404      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);
405      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);
406      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);
407      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);
408      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);
409      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);
410      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);
411      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);
412      //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);
413      //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);
414      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);
415      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);
416      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);
417      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);
418      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);
419      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);
420      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);
421      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);
422      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);
423      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);
424      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);
425      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);
426      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);
427      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);
428      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);
429      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);
430      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);
431      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);
432      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);
433      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);
434      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);
435      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);
436      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);
437      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);
438      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);
439      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);
440      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);
441      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);
442      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);
443      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);
444      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);
445      //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);
446      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);
447      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);
448      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);
449      //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);
450      //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);
451      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);
452      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);
453      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);
454      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);
455    } else {
456      // Just testing one
457      mpsName.push_back(empty.problemName());min.push_back(true);nRows.push_back(-1);
458      nCols.push_back(-1);objValueTol.push_back(1.e-10);
459      objValue.push_back(0.0);bestStrategy.push_back(0);
460      int iTest;
461      std::string alg;
462      for (iTest=0;iTest<NUMBER_ALGORITHMS;iTest++) {
463        ClpSolve solveOptions=setupForSolve(iTest,alg,0);
464        printf("%d %s ",iTest,alg.c_str());
465        if (switchOff[iTest]) 
466          printf("skipped by user\n");
467        else if(solveOptions.getSolveType()==ClpSolve::notImplemented)
468          printf("skipped as not available\n");
469        else
470          printf("will be tested\n");
471      }
472    }
473
474    double timeTaken =0.0;
475    if( !barrierAvailable)
476      switchOff[0]=1;
477  // Loop once for each Mps File
478    for (m=0; m<mpsName.size(); m++ ) {
479      std::cerr <<"  processing mps file: " <<mpsName[m] 
480                <<" (" <<m+1 <<" out of " <<mpsName.size() <<")" <<std::endl;
481
482      ClpSimplex solutionBase=empty;
483      std::string fn = dirNetlib+mpsName[m];
484      if (!empty.numberRows()||algorithm<6) {
485        // Read data mps file,
486        CoinMpsIO mps;
487        int nerrors=mps.readMps(fn.c_str(),"mps");
488        if (nerrors) {
489          std::cerr << "Error " << nerrors << " when reading model from "
490                    << fn.c_str() << "! Aborting tests.\n";
491          return 1;
492        }
493        solutionBase.loadProblem(*mps.getMatrixByCol(),mps.getColLower(),
494                                 mps.getColUpper(),
495                                 mps.getObjCoefficients(),
496                                 mps.getRowLower(),mps.getRowUpper());
497       
498        solutionBase.setDblParam(ClpObjOffset,mps.objectiveOffset());
499      } 
500     
501      // Runs through strategies
502      if (algorithm==6||algorithm==7) {
503        // algorithms tested are at top of file
504        double testTime[NUMBER_ALGORITHMS];
505        std::string alg[NUMBER_ALGORITHMS];
506        int iTest;
507        for (iTest=0;iTest<NUMBER_ALGORITHMS;iTest++) {
508          ClpSolve solveOptions=setupForSolve(iTest,alg[iTest],1);
509          if (solveOptions.getSolveType()!=ClpSolve::notImplemented) {
510            double time1 = CoinCpuTime();
511            ClpSimplex solution=solutionBase;
512            if (solution.maximumSeconds()<0.0)
513              solution.setMaximumSeconds(120.0);
514            if (doVector) {
515              ClpMatrixBase * matrix = solution.clpMatrix();
516              if (dynamic_cast< ClpPackedMatrix*>(matrix)) {
517                ClpPackedMatrix * clpMatrix = dynamic_cast< ClpPackedMatrix*>(matrix);
518                clpMatrix->makeSpecialColumnCopy();
519              }
520            }
521            solution.initialSolve(solveOptions);
522            double time2 = CoinCpuTime()-time1;
523            testTime[iTest]=time2;
524            printf("Took %g seconds - status %d\n",time2,solution.problemStatus());
525            if (solution.problemStatus()) 
526              testTime[iTest]=1.0e20;
527          } else {
528            testTime[iTest]=1.0e30;
529          }
530        }
531        int iBest=-1;
532        double dBest=1.0e10;
533        printf("%s",fn.c_str());
534        for (iTest=0;iTest<NUMBER_ALGORITHMS;iTest++) {
535          if (testTime[iTest]<1.0e30) {
536            printf(" %s %g",
537                   alg[iTest].c_str(),testTime[iTest]);
538            if (testTime[iTest]<dBest) {
539              dBest=testTime[iTest];
540              iBest=iTest;
541            }
542          }
543        }
544        printf("\n");
545        if (iBest>=0)
546          printf("Best strategy for %s is %s (%d) which takes %g seconds\n",
547                 fn.c_str(),alg[iBest].c_str(),iBest,testTime[iBest]);
548        else
549          printf("No strategy finished in time limit\n");
550        continue;
551      }
552      double time1 = CoinCpuTime();
553      ClpSimplex solution=solutionBase;
554#if 0
555      solution.setOptimizationDirection(-1);
556      {
557        int j;
558        double * obj = solution.objective();
559        int n=solution.numberColumns();
560        for (j=0;j<n;j++)
561          obj[j] *= -1.0;
562      }
563#endif
564      ClpSolve::SolveType method;
565      ClpSolve solveOptions=solveOptionsIn;
566      std::string nameAlgorithm;
567      if (algorithm!=5) {
568        if (algorithm==0) {
569          method=ClpSolve::useDual;
570          nameAlgorithm="dual";
571        } else if (algorithm==1) {
572          method=ClpSolve::usePrimalorSprint;
573          nameAlgorithm="primal";
574        } else if (algorithm==3) {
575          method=ClpSolve::automatic;
576          nameAlgorithm="either";
577        } else {
578          nameAlgorithm="barrier-slow";
579#ifdef UFL_BARRIER
580          solveOptions.setSpecialOption(4,4);
581          nameAlgorithm="barrier-UFL";
582#endif
583#ifdef WSSMP_BARRIER
584          solveOptions.setSpecialOption(4,2);
585          nameAlgorithm="barrier-WSSMP";
586#endif
587#ifdef MUMPS_BARRIER
588          solveOptions.setSpecialOption(4,6);
589          nameAlgorithm="barrier-MUMPS";
590#endif
591          method = ClpSolve::useBarrier;
592        }
593        solveOptions.setSolveType(method);
594      } else {
595        int iAlg = bestStrategy[m];
596        int presolveOff=iAlg/100;
597        iAlg=iAlg % 100;
598        if( !barrierAvailable&&iAlg==0) {
599          if (nRows[m]!=2172)
600            iAlg = 5; // try primal
601          else
602            iAlg=3; // d2q06c
603        }
604        solveOptions=setupForSolve(iAlg,nameAlgorithm,0);
605        if (presolveOff)
606          solveOptions.setPresolveType(ClpSolve::presolveOff);
607      }
608      if (doVector) {
609        ClpMatrixBase * matrix = solution.clpMatrix();
610        if (dynamic_cast< ClpPackedMatrix*>(matrix)) {
611          ClpPackedMatrix * clpMatrix = dynamic_cast< ClpPackedMatrix*>(matrix);
612          clpMatrix->makeSpecialColumnCopy();
613        }
614      }
615      solution.initialSolve(solveOptions);
616      double time2 = CoinCpuTime()-time1;
617      timeTaken += time2;
618      printf("%s took %g seconds using algorithm %s\n",fn.c_str(),time2,nameAlgorithm.c_str());
619      // Test objective solution value
620      {
621        double soln = solution.objectiveValue();
622        CoinRelFltEq eq(objValueTol[m]);
623        std::cerr <<soln <<",  " <<objValue[m] <<" diff "<<
624          soln-objValue[m]<<std::endl;
625        if(!eq(soln,objValue[m]))
626          printf("** difference fails\n");
627      }
628    }
629    printf("Total time %g seconds\n",timeTaken);
630  }
631  else {
632    testingMessage( "***Skipped Testing on netlib    ***\n" );
633    testingMessage( "***use -netlib to test class***\n" );
634  }
635 
636  testingMessage( "All tests completed successfully\n" );
637  return 0;
638}
639
640 
641// Display message on stdout and stderr
642void testingMessage( const char * const msg )
643{
644  std::cerr <<msg;
645  //cout <<endl <<"*****************************************"
646  //     <<endl <<msg <<endl;
647}
648
649//--------------------------------------------------------------------------
650// test factorization methods and simplex method and simple barrier
651void
652ClpSimplexUnitTest(const std::string & dirSample)
653{
654 
655  CoinRelFltEq eq(0.000001);
656
657  {
658    ClpSimplex solution;
659 
660    // matrix data
661    //deliberate hiccup of 2 between 0 and 1
662    CoinBigIndex start[5]={0,4,7,8,9};
663    int length[5]={2,3,1,1,1};
664    int rows[11]={0,2,-1,-1,0,1,2,0,1,2};
665    double elements[11]={7.0,2.0,1.0e10,1.0e10,-2.0,1.0,-2.0,1,1,1};
666    CoinPackedMatrix matrix(true,3,5,8,elements,rows,start,length);
667   
668    // rim data
669    double objective[7]={-4.0,1.0,0.0,0.0,0.0,0.0,0.0};
670    double rowLower[5]={14.0,3.0,3.0,1.0e10,1.0e10};
671    double rowUpper[5]={14.0,3.0,3.0,-1.0e10,-1.0e10};
672    double colLower[7]={0.0,0.0,0.0,0.0,0.0,0.0,0.0};
673    double colUpper[7]={100.0,100.0,100.0,100.0,100.0,100.0,100.0};
674   
675    // basis 1
676    int rowBasis1[3]={-1,-1,-1};
677    int colBasis1[5]={1,1,-1,-1,1};
678    solution.loadProblem(matrix,colLower,colUpper,objective,
679                         rowLower,rowUpper);
680    int i;
681    solution.createStatus();
682    for (i=0;i<3;i++) {
683      if (rowBasis1[i]<0) {
684        solution.setRowStatus(i,ClpSimplex::atLowerBound);
685      } else {
686        solution.setRowStatus(i,ClpSimplex::basic);
687      }
688    }
689    for (i=0;i<5;i++) {
690      if (colBasis1[i]<0) {
691        solution.setColumnStatus(i,ClpSimplex::atLowerBound);
692      } else {
693        solution.setColumnStatus(i,ClpSimplex::basic);
694      }
695    }
696    solution.setLogLevel(3+4+8+16+32);
697    solution.primal();
698    for (i=0;i<3;i++) {
699      if (rowBasis1[i]<0) {
700        solution.setRowStatus(i,ClpSimplex::atLowerBound);
701      } else {
702        solution.setRowStatus(i,ClpSimplex::basic);
703      }
704    }
705    for (i=0;i<5;i++) {
706      if (colBasis1[i]<0) {
707        solution.setColumnStatus(i,ClpSimplex::atLowerBound);
708      } else {
709        solution.setColumnStatus(i,ClpSimplex::basic);
710      }
711    }
712    // intricate stuff does not work with scaling
713    solution.scaling(0);
714#ifndef NDEBUG
715    int returnCode = solution.factorize ( );
716    assert(!returnCode);
717#else
718    solution.factorize ( );
719#endif
720    const double * colsol = solution.primalColumnSolution();
721    const double * rowsol = solution.primalRowSolution();
722    solution.getSolution(rowsol,colsol);
723#ifndef NDEBUG
724    double colsol1[5]={20.0/7.0,3.0,0.0,0.0,23.0/7.0};
725    for (i=0;i<5;i++) {
726      assert(eq(colsol[i],colsol1[i]));
727    }
728    // now feed in again without actually doing factorization
729#endif
730    ClpFactorization factorization2 = *solution.factorization();
731    ClpSimplex solution2 = solution;
732    solution2.setFactorization(factorization2);
733    solution2.createStatus();
734    for (i=0;i<3;i++) {
735      if (rowBasis1[i]<0) {
736        solution2.setRowStatus(i,ClpSimplex::atLowerBound);
737      } else {
738        solution2.setRowStatus(i,ClpSimplex::basic);
739      }
740    }
741    for (i=0;i<5;i++) {
742      if (colBasis1[i]<0) {
743        solution2.setColumnStatus(i,ClpSimplex::atLowerBound);
744      } else {
745        solution2.setColumnStatus(i,ClpSimplex::basic);
746      }
747    }
748    // intricate stuff does not work with scaling
749    solution2.scaling(0);
750    solution2.getSolution(rowsol,colsol);
751    colsol = solution2.primalColumnSolution();
752    rowsol = solution2.primalRowSolution();
753    for (i=0;i<5;i++) {
754      assert(eq(colsol[i],colsol1[i]));
755    }
756    solution2.setDualBound(0.1);
757    solution2.dual();
758    objective[2]=-1.0;
759    objective[3]=-0.5;
760    objective[4]=10.0;
761    solution.dual();
762    for (i=0;i<3;i++) {
763      rowLower[i]=-1.0e20;
764      colUpper[i+2]=0.0;
765    }
766    solution.setLogLevel(3);
767    solution.dual();
768    double rowObjective[]={1.0,0.5,-10.0};
769    solution.loadProblem(matrix,colLower,colUpper,objective,
770                         rowLower,rowUpper,rowObjective);
771    solution.dual();
772    solution.loadProblem(matrix,colLower,colUpper,objective,
773                         rowLower,rowUpper,rowObjective);
774    solution.primal();
775  }
776#ifndef COIN_NO_CLP_MESSAGE
777  {   
778    CoinMpsIO m;
779    std::string fn = dirSample+"exmip1";
780    m.readMps(fn.c_str(),"mps");
781    ClpSimplex solution;
782    solution.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
783                         m.getObjCoefficients(),
784                         m.getRowLower(),m.getRowUpper());
785    solution.dual();
786    // Test event handling
787    MyEventHandler handler;
788    solution.passInEventHandler(&handler);
789    int numberRows=solution.numberRows();
790    // make sure values pass has something to do
791    for (int i=0;i<numberRows;i++)
792      solution.setRowStatus(i,ClpSimplex::basic);
793    solution.primal(1);
794    assert (solution.secondaryStatus()==102); // Came out at end of pass
795  }
796  // Test Message handler
797  {   
798    CoinMpsIO m;
799    std::string fn = dirSample+"exmip1";
800    //fn = "Test/subGams4";
801    m.readMps(fn.c_str(),"mps");
802    ClpSimplex model;
803    model.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
804                         m.getObjCoefficients(),
805                         m.getRowLower(),m.getRowUpper());
806    // Message handler
807    MyMessageHandler messageHandler(&model);
808    std::cout<<"Testing derived message handler"<<std::endl;
809    model.passInMessageHandler(&messageHandler);
810    model.messagesPointer()->setDetailMessage(1,102);
811    model.setFactorizationFrequency(10);
812    model.primal();
813    model.primal(0,3);
814    model.setObjCoeff(3,-2.9473684210526314);
815    model.primal(0,3);
816    // Write saved solutions
817    int nc = model.getNumCols();
818    int s; 
819    std::deque<StdVectorDouble> fep = messageHandler.getFeasibleExtremePoints();
820    int numSavedSolutions = fep.size();
821    for ( s=0; s<numSavedSolutions; ++s ) {
822      const StdVectorDouble & solnVec = fep[s];
823      for ( int c=0; c<nc; ++c ) {
824        if (fabs(solnVec[c])>1.0e-8)
825          std::cout <<"Saved Solution: " <<s <<" ColNum: " <<c <<" Value: " <<solnVec[c] <<std::endl;
826      }
827    }
828    // Solve again without scaling
829    // and maximize then minimize
830    messageHandler.clearFeasibleExtremePoints();
831    model.scaling(0);
832    model.setOptimizationDirection(-1);
833    model.primal();
834    model.setOptimizationDirection(1);
835    model.primal();
836    fep = messageHandler.getFeasibleExtremePoints();
837    numSavedSolutions = fep.size();
838    for ( s=0; s<numSavedSolutions; ++s ) {
839      const StdVectorDouble & solnVec = fep[s];
840      for ( int c=0; c<nc; ++c ) {
841        if (fabs(solnVec[c])>1.0e-8)
842          std::cout <<"Saved Solution: " <<s <<" ColNum: " <<c <<" Value: " <<solnVec[c] <<std::endl;
843      }
844    }
845  }
846#endif
847  // Test dual ranging
848  {   
849    CoinMpsIO m;
850    std::string fn = dirSample+"exmip1";
851    m.readMps(fn.c_str(),"mps");
852    ClpSimplex model;
853    model.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
854                         m.getObjCoefficients(),
855                         m.getRowLower(),m.getRowUpper());
856    model.primal();
857    int which[13] = {0,1,2,3,4,5,6,7,8,9,10,11,12};
858    double costIncrease[13];
859    int sequenceIncrease[13];
860    double costDecrease[13];
861    int sequenceDecrease[13];
862    // ranging
863    model.dualRanging(13,which,costIncrease,sequenceIncrease,
864                      costDecrease,sequenceDecrease);
865    int i;
866    for ( i=0;i<13;i++)
867      printf("%d increase %g %d, decrease %g %d\n",
868             i,costIncrease[i],sequenceIncrease[i],
869             costDecrease[i],sequenceDecrease[i]);
870    assert (fabs(costDecrease[3])<1.0e-4);
871    assert (fabs(costIncrease[7]-1.0)<1.0e-4);
872    model.setOptimizationDirection(-1);
873    {
874      int j;
875      double * obj = model.objective();
876      int n=model.numberColumns();
877      for (j=0;j<n;j++) 
878        obj[j] *= -1.0;
879    }
880    double costIncrease2[13];
881    int sequenceIncrease2[13];
882    double costDecrease2[13];
883    int sequenceDecrease2[13];
884    // ranging
885    model.dualRanging(13,which,costIncrease2,sequenceIncrease2,
886                      costDecrease2,sequenceDecrease2);
887    for (i=0;i<13;i++) {
888      assert (fabs(costIncrease[i]-costDecrease2[i])<1.0e-6);
889      assert (fabs(costDecrease[i]-costIncrease2[i])<1.0e-6);
890      assert (sequenceIncrease[i]==sequenceDecrease2[i]);
891      assert (sequenceDecrease[i]==sequenceIncrease2[i]);
892    }
893    // Now delete all rows and see what happens
894    model.deleteRows(model.numberRows(),which);
895    model.primal();
896    // ranging
897    if (!model.dualRanging(8,which,costIncrease,sequenceIncrease,
898                           costDecrease,sequenceDecrease)) {
899      for (i=0;i<8;i++) {
900        printf("%d increase %g %d, decrease %g %d\n",
901               i,costIncrease[i],sequenceIncrease[i],
902               costDecrease[i],sequenceDecrease[i]);
903      }
904    }
905  }
906  // Test primal ranging
907  {   
908    CoinMpsIO m;
909    std::string fn = dirSample+"exmip1";
910    m.readMps(fn.c_str(),"mps");
911    ClpSimplex model;
912    model.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
913                         m.getObjCoefficients(),
914                         m.getRowLower(),m.getRowUpper());
915    model.primal();
916    int which[13] = {0,1,2,3,4,5,6,7,8,9,10,11,12};
917    double valueIncrease[13];
918    int sequenceIncrease[13];
919    double valueDecrease[13];
920    int sequenceDecrease[13];
921    // ranging
922    model.primalRanging(13,which,valueIncrease,sequenceIncrease,
923                      valueDecrease,sequenceDecrease);
924    int i;
925    for ( i=0;i<13;i++)
926      printf("%d increase %g %d, decrease %g %d\n",
927             i,valueIncrease[i],sequenceIncrease[i],
928             valueDecrease[i],sequenceDecrease[i]);
929    assert (fabs(valueIncrease[3]-0.642857)<1.0e-4);
930    assert (fabs(valueIncrease[8]-2.95113)<1.0e-4);
931#if 0
932    // out until I find optimization bug
933    // Test parametrics
934    ClpSimplexOther * model2 = (ClpSimplexOther *) (&model);
935    double rhs[]={ 1.0,2.0,3.0,4.0,5.0};
936    double endingTheta=1.0;
937    model2->scaling(0);
938    model2->setLogLevel(63);
939    model2->parametrics(0.0,endingTheta,0.1,
940                        NULL,NULL,rhs,rhs,NULL);
941#endif
942  }
943  // Test binv etc
944  {   
945    /*
946       Wolsey : Page 130
947       max 4x1 -  x2
948       7x1 - 2x2    <= 14
949       x2    <= 3
950       2x1 - 2x2    <= 3
951       x1 in Z+, x2 >= 0
952
953       note slacks are -1 in Clp so signs may be different
954    */
955   
956    int n_cols = 2;
957    int n_rows = 3;
958   
959    double obj[2] = {-4.0, 1.0};
960    double collb[2] = {0.0, 0.0};
961    double colub[2] = {COIN_DBL_MAX, COIN_DBL_MAX};
962    double rowlb[3] = {-COIN_DBL_MAX, -COIN_DBL_MAX, -COIN_DBL_MAX};
963    double rowub[3] = {14.0, 3.0, 3.0};
964   
965    int rowIndices[5] =  {0,     2,    0,    1,    2};
966    int colIndices[5] =  {0,     0,    1,    1,    1};
967    double elements[5] = {7.0, 2.0, -2.0,  1.0, -2.0};
968    CoinPackedMatrix M(true, rowIndices, colIndices, elements, 5);
969
970    ClpSimplex model;
971    model.loadProblem(M, collb, colub, obj, rowlb, rowub);
972    model.dual(0,1); // keep factorization
973   
974    //check that the tableau matches wolsey (B-1 A)
975    // slacks in second part of binvA
976    double * binvA = reinterpret_cast<double*> (malloc((n_cols+n_rows) * sizeof(double)));
977   
978    printf("B-1 A by row\n");
979    int i;
980    for( i = 0; i < n_rows; i++){
981      model.getBInvARow(i, binvA,binvA+n_cols);
982      printf("row: %d -> ",i);
983      for(int j=0; j < n_cols+n_rows; j++){
984        printf("%g, ", binvA[j]);
985      }
986      printf("\n");
987    }
988    // See if can re-use factorization AND arrays
989    model.primal(0,3+4); // keep factorization
990    // And do by column
991    printf("B-1 A by column\n");
992    for( i = 0; i < n_rows+n_cols; i++){
993      model.getBInvACol(i, binvA);
994      printf("column: %d -> ",i);
995      for(int j=0; j < n_rows; j++){
996        printf("%g, ", binvA[j]);
997      }
998      printf("\n");
999    }
1000    /* Do twice -
1001       without and with scaling
1002    */
1003    // set scaling off
1004    model.scaling(0);
1005    for (int iPass=0;iPass<2;iPass++) {
1006      model.primal(0,3+4); // keep factorization
1007      const double * rowScale = model.rowScale();
1008      const double * columnScale = model.columnScale();
1009      if (!iPass)
1010        assert (!rowScale);
1011      else
1012        assert (rowScale); // only true for this example
1013      /* has to be exactly correct as in OsiClpsolverInterface.cpp
1014         (also redo each pass as may change
1015      */
1016      printf("B-1 A");
1017      for( i = 0; i < n_rows; i++){
1018        model.getBInvARow(i, binvA,binvA+n_cols);
1019        printf("\nrow: %d -> ",i);
1020        int j;
1021        // First columns
1022        for(j=0; j < n_cols; j++){
1023          if (binvA[j]) {
1024            printf("(%d %g), ", j, binvA[j]);
1025          }
1026        }
1027        // now rows
1028        for(j=0; j < n_rows; j++){
1029          if (binvA[j+n_cols]) {
1030            printf("(%d %g), ", j+n_cols, binvA[j+n_cols]);
1031          }
1032        }
1033      }
1034      printf("\n");
1035      printf("And by column (trickier)");
1036      const int * pivotVariable = model.pivotVariable();
1037      for( i = 0; i < n_cols+n_rows; i++){
1038        model.getBInvACol(i, binvA);
1039        printf("\ncolumn: %d -> ",i);
1040        for(int j=0; j < n_rows; j++){
1041          if (binvA[j]) {
1042            // need to know pivot variable for +1/-1 (slack) and row/column scaling
1043            int pivot = pivotVariable[j];
1044            if (pivot<n_cols) {
1045              // scaled coding is in just in case
1046              if (!columnScale) {
1047                printf("(%d %g), ", j, binvA[j]);
1048              } else {
1049                printf("(%d %g), ", j, binvA[j]*columnScale[pivot]);
1050              }
1051            } else {
1052              if (!rowScale) {
1053                printf("(%d %g), ", j, binvA[j]);
1054              } else {
1055                printf("(%d %g), ", j, binvA[j]/rowScale[pivot-n_cols]);
1056              }
1057            }
1058          }
1059        }
1060      }
1061      printf("\n");
1062      printf("binvrow");
1063      for( i = 0; i < n_rows; i++){
1064        model.getBInvRow(i, binvA);
1065        printf("\nrow: %d -> ",i);
1066        int j;
1067        for (j=0;j<n_rows;j++) {
1068          if (binvA[j])
1069            printf("(%d %g), ", j, binvA[j]);
1070        }
1071      }
1072      printf("\n");
1073      printf("And by column ");
1074      for( i = 0; i < n_rows; i++){
1075        model.getBInvCol(i, binvA);
1076        printf("\ncol: %d -> ",i);
1077        int j;
1078        for (j=0;j<n_rows;j++) {
1079          if (binvA[j])
1080            printf("(%d %g), ", j, binvA[j]);
1081        }
1082      }
1083      printf("\n");
1084      // now deal with next pass
1085      if (!iPass) {
1086        // get scaling for testing
1087        model.scaling(1);
1088      }
1089    }
1090    free(binvA);
1091    model.setColUpper(1,2.0);
1092    model.dual(0,2+4); // use factorization and arrays
1093    model.dual(0,2); // hopefully will not use factorization
1094    model.primal(0,3+4); // keep factorization
1095    // but say basis has changed
1096    model.setWhatsChanged(model.whatsChanged()&(~512));
1097    model.dual(0,2); // hopefully will not use factorization
1098  }
1099  // test steepest edge
1100  {   
1101    CoinMpsIO m;
1102    std::string fn = dirSample+"finnis";
1103    int returnCode = m.readMps(fn.c_str(),"mps");
1104    if (returnCode) {
1105      // probable cause is that gz not there
1106      fprintf(stderr,"Unable to open finnis.mps in %s\n", dirSample.c_str());
1107      fprintf(stderr,"Most probable cause is finnis.mps is gzipped i.e. finnis.mps.gz and libz has not been activated\n");
1108      fprintf(stderr,"Either gunzip files or edit Makefiles/Makefile.location to get libz\n");
1109      exit(999);
1110    }
1111    ClpModel model;
1112    model.loadProblem(*m.getMatrixByCol(),m.getColLower(),
1113                    m.getColUpper(),
1114                    m.getObjCoefficients(),
1115                    m.getRowLower(),m.getRowUpper());
1116    ClpSimplex solution(model);
1117
1118    solution.scaling(1); 
1119    solution.setDualBound(1.0e8);
1120    //solution.factorization()->maximumPivots(1);
1121    //solution.setLogLevel(3);
1122    solution.setDualTolerance(1.0e-7);
1123    // set objective sense,
1124    ClpDualRowSteepest steep;
1125    solution.setDualRowPivotAlgorithm(steep);
1126    solution.setDblParam(ClpObjOffset,m.objectiveOffset());
1127    solution.dual();
1128  }
1129  // test normal solution
1130  {   
1131    CoinMpsIO m;
1132    std::string fn = dirSample+"afiro";
1133    m.readMps(fn.c_str(),"mps");
1134    ClpSimplex solution;
1135    ClpModel model;
1136    // do twice - without and with scaling
1137    int iPass;
1138    for (iPass=0;iPass<2;iPass++) {
1139      // explicit row objective for testing
1140      int nr = m.getNumRows();
1141      double * rowObj = new double[nr];
1142      CoinFillN(rowObj,nr,0.0);
1143      model.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
1144                      m.getObjCoefficients(),
1145                      m.getRowLower(),m.getRowUpper(),rowObj);
1146      delete [] rowObj;
1147      solution = ClpSimplex(model);
1148      if (iPass) {
1149        solution.scaling();
1150      }
1151      solution.dual();
1152      solution.dual();
1153      // test optimal
1154      assert (solution.status()==0);
1155      int numberColumns = solution.numberColumns();
1156      int numberRows = solution.numberRows();
1157      CoinPackedVector colsol(numberColumns,solution.primalColumnSolution());
1158      double * objective = solution.objective();
1159#ifndef NDEBUG
1160      double objValue = colsol.dotProduct(objective);
1161#endif
1162      CoinRelFltEq eq(1.0e-8);
1163      assert(eq(objValue,-4.6475314286e+02));
1164      solution.dual();
1165      assert(eq(solution.objectiveValue(),-4.6475314286e+02));
1166      double * lower = solution.columnLower();
1167      double * upper = solution.columnUpper();
1168      double * sol = solution.primalColumnSolution();
1169      double * result = new double[numberColumns];
1170      CoinFillN ( result, numberColumns,0.0);
1171      solution.matrix()->transposeTimes(solution.dualRowSolution(), result);
1172      int iRow , iColumn;
1173      // see if feasible and dual feasible
1174      for (iColumn=0;iColumn<numberColumns;iColumn++) {
1175        double value = sol[iColumn];
1176        assert(value<upper[iColumn]+1.0e-8);
1177        assert(value>lower[iColumn]-1.0e-8);
1178        value = objective[iColumn]-result[iColumn];
1179        assert (value>-1.0e-5);
1180        if (sol[iColumn]>1.0e-5)
1181          assert (value<1.0e-5);
1182      }
1183      delete [] result;
1184      result = new double[numberRows];
1185      CoinFillN ( result, numberRows,0.0);
1186      solution.matrix()->times(colsol, result);
1187      lower = solution.rowLower();
1188      upper = solution.rowUpper();
1189      sol = solution.primalRowSolution();
1190#ifndef NDEBUG
1191      for (iRow=0;iRow<numberRows;iRow++) {
1192        double value = result[iRow];
1193        assert(eq(value,sol[iRow]));
1194        assert(value<upper[iRow]+1.0e-8);
1195        assert(value>lower[iRow]-1.0e-8);
1196      }
1197#endif
1198      delete [] result;
1199      // test row objective
1200      double * rowObjective = solution.rowObjective();
1201      CoinDisjointCopyN(solution.dualRowSolution(),numberRows,rowObjective);
1202      CoinDisjointCopyN(solution.dualColumnSolution(),numberColumns,objective);
1203      // this sets up all slack basis
1204      solution.createStatus();
1205      solution.dual();
1206      CoinFillN(rowObjective,numberRows,0.0);
1207      CoinDisjointCopyN(m.getObjCoefficients(),numberColumns,objective);
1208      solution.dual();
1209    }
1210  }
1211  // test unbounded
1212  {   
1213    CoinMpsIO m;
1214    std::string fn = dirSample+"brandy";
1215    m.readMps(fn.c_str(),"mps");
1216    ClpSimplex solution;
1217    // do twice - without and with scaling
1218    int iPass;
1219    for (iPass=0;iPass<2;iPass++) {
1220      solution.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
1221                      m.getObjCoefficients(),
1222                      m.getRowLower(),m.getRowUpper());
1223      if (iPass)
1224        solution.scaling();
1225      solution.setOptimizationDirection(-1);
1226      // test unbounded and ray
1227#ifdef DUAL
1228      solution.setDualBound(100.0);
1229      solution.dual();
1230#else
1231      solution.primal();
1232#endif
1233      assert (solution.status()==2);
1234      int numberColumns = solution.numberColumns();
1235      int numberRows = solution.numberRows();
1236      double * lower = solution.columnLower();
1237      double * upper = solution.columnUpper();
1238      double * sol = solution.primalColumnSolution();
1239      double * ray = solution.unboundedRay();
1240      double * objective = solution.objective();
1241      double objChange=0.0;
1242      int iRow , iColumn;
1243      // make sure feasible and columns form ray
1244      for (iColumn=0;iColumn<numberColumns;iColumn++) {
1245        double value = sol[iColumn];
1246        assert(value<upper[iColumn]+1.0e-8);
1247        assert(value>lower[iColumn]-1.0e-8);
1248        value = ray[iColumn];
1249        if (value>0.0)
1250          assert(upper[iColumn]>1.0e30);
1251        else if (value<0.0)
1252          assert(lower[iColumn]<-1.0e30);
1253        objChange += value*objective[iColumn];
1254      }
1255      // make sure increasing objective
1256      assert(objChange>0.0);
1257      double * result = new double[numberRows];
1258      CoinFillN ( result, numberRows,0.0);
1259      solution.matrix()->times(sol, result);
1260      lower = solution.rowLower();
1261      upper = solution.rowUpper();
1262      sol = solution.primalRowSolution();
1263#ifndef NDEBUG
1264      for (iRow=0;iRow<numberRows;iRow++) {
1265        double value = result[iRow];
1266        assert(eq(value,sol[iRow]));
1267        assert(value<upper[iRow]+2.0e-8);
1268        assert(value>lower[iRow]-2.0e-8);
1269      }
1270#endif
1271      CoinFillN ( result, numberRows,0.0);
1272      solution.matrix()->times(ray, result);
1273      // there may be small differences (especially if scaled)
1274      for (iRow=0;iRow<numberRows;iRow++) {
1275        double value = result[iRow];
1276        if (value>1.0e-8)
1277          assert(upper[iRow]>1.0e30);
1278        else if (value<-1.0e-8)
1279          assert(lower[iRow]<-1.0e30);
1280      }
1281      delete [] result;
1282      delete [] ray;
1283    }
1284  }
1285  // test infeasible
1286  {   
1287    CoinMpsIO m;
1288    std::string fn = dirSample+"brandy";
1289    m.readMps(fn.c_str(),"mps");
1290    ClpSimplex solution;
1291    // do twice - without and with scaling
1292    int iPass;
1293    for (iPass=0;iPass<2;iPass++) {
1294      solution.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
1295                      m.getObjCoefficients(),
1296                      m.getRowLower(),m.getRowUpper());
1297      if (iPass)
1298        solution.scaling();
1299      // test infeasible and ray
1300      solution.columnUpper()[0]=0.0;
1301#ifdef DUAL
1302      solution.setDualBound(100.0);
1303      solution.dual();
1304#else
1305      solution.primal();
1306#endif
1307      assert (solution.status()==1);
1308      int numberColumns = solution.numberColumns();
1309      int numberRows = solution.numberRows();
1310      double * lower = solution.rowLower();
1311      double * upper = solution.rowUpper();
1312      double * ray = solution.infeasibilityRay();
1313      assert(ray);
1314      // construct proof of infeasibility
1315      int iRow , iColumn;
1316      double lo=0.0,up=0.0;
1317      int nl=0,nu=0;
1318      for (iRow=0;iRow<numberRows;iRow++) {
1319        if (lower[iRow]>-1.0e20) {
1320          lo += ray[iRow]*lower[iRow];
1321        } else {
1322          if (ray[iRow]>1.0e-8) 
1323            nl++;
1324        }
1325        if (upper[iRow]<1.0e20) {
1326          up += ray[iRow]*upper[iRow];
1327        } else {
1328          if (ray[iRow]>1.0e-8) 
1329            nu++;
1330        }
1331      }
1332      if (nl)
1333        lo=-1.0e100;
1334      if (nu)
1335        up=1.0e100;
1336      double * result = new double[numberColumns];
1337      double lo2=0.0,up2=0.0;
1338      CoinFillN ( result, numberColumns,0.0);
1339      solution.matrix()->transposeTimes(ray, result);
1340      lower = solution.columnLower();
1341      upper = solution.columnUpper();
1342      nl=nu=0;
1343      for (iColumn=0;iColumn<numberColumns;iColumn++) {
1344        if (result[iColumn]>1.0e-8) {
1345          if (lower[iColumn]>-1.0e20)
1346            lo2 += result[iColumn]*lower[iColumn];
1347          else
1348            nl++;
1349          if (upper[iColumn]<1.0e20)
1350            up2 += result[iColumn]*upper[iColumn];
1351          else
1352            nu++;
1353        } else if (result[iColumn]<-1.0e-8) {
1354          if (lower[iColumn]>-1.0e20)
1355            up2 += result[iColumn]*lower[iColumn];
1356          else
1357            nu++;
1358          if (upper[iColumn]<1.0e20)
1359            lo2 += result[iColumn]*upper[iColumn];
1360          else
1361            nl++;
1362        }
1363      }
1364      if (nl)
1365        lo2=-1.0e100;
1366      if (nu)
1367        up2=1.0e100;
1368      // make sure inconsistency
1369      assert(lo2>up||up2<lo);
1370      delete [] result;
1371      delete [] ray;
1372    }
1373  }
1374  // test delete and add
1375  {   
1376    CoinMpsIO m;
1377    std::string fn = dirSample+"brandy";
1378    m.readMps(fn.c_str(),"mps");
1379    ClpSimplex solution;
1380    solution.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
1381                         m.getObjCoefficients(),
1382                         m.getRowLower(),m.getRowUpper());
1383    solution.dual();
1384    CoinRelFltEq eq(1.0e-8);
1385    assert(eq(solution.objectiveValue(),1.5185098965e+03));
1386
1387    int numberColumns = solution.numberColumns();
1388    int numberRows = solution.numberRows();
1389    double * saveObj = new double [numberColumns];
1390    double * saveLower = new double[numberRows+numberColumns];
1391    double * saveUpper = new double[numberRows+numberColumns];
1392    int * which = new int [numberRows+numberColumns];
1393
1394    int numberElements = m.getMatrixByCol()->getNumElements();
1395    int * starts = new int[numberRows+numberColumns];
1396    int * index = new int[numberElements];
1397    double * element = new double[numberElements];
1398
1399    const CoinBigIndex * startM;
1400    const int * lengthM;
1401    const int * indexM;
1402    const double * elementM;
1403
1404    int n,nel;
1405
1406    // delete non basic columns
1407    n=0;
1408    nel=0;
1409    int iRow , iColumn;
1410    const double * lower = m.getColLower();
1411    const double * upper = m.getColUpper();
1412    const double * objective = m.getObjCoefficients();
1413    startM = m.getMatrixByCol()->getVectorStarts();
1414    lengthM = m.getMatrixByCol()->getVectorLengths();
1415    indexM = m.getMatrixByCol()->getIndices();
1416    elementM = m.getMatrixByCol()->getElements();
1417    starts[0]=0;
1418    for (iColumn=0;iColumn<numberColumns;iColumn++) {
1419      if (solution.getColumnStatus(iColumn)!=ClpSimplex::basic) {
1420        saveObj[n]=objective[iColumn];
1421        saveLower[n]=lower[iColumn];
1422        saveUpper[n]=upper[iColumn];
1423        int j;
1424        for (j=startM[iColumn];j<startM[iColumn]+lengthM[iColumn];j++) {
1425          index[nel]=indexM[j];
1426          element[nel++]=elementM[j];
1427        }
1428        which[n++]=iColumn;
1429        starts[n]=nel;
1430      }
1431    }
1432    solution.deleteColumns(n,which);
1433    solution.dual();
1434    // Put back
1435    solution.addColumns(n,saveLower,saveUpper,saveObj,
1436                        starts,index,element);
1437    solution.dual();
1438    assert(eq(solution.objectiveValue(),1.5185098965e+03));
1439    // Delete all columns and add back
1440    n=0;
1441    nel=0;
1442    starts[0]=0;
1443    lower = m.getColLower();
1444    upper = m.getColUpper();
1445    objective = m.getObjCoefficients();
1446    for (iColumn=0;iColumn<numberColumns;iColumn++) {
1447      saveObj[n]=objective[iColumn];
1448      saveLower[n]=lower[iColumn];
1449      saveUpper[n]=upper[iColumn];
1450      int j;
1451      for (j=startM[iColumn];j<startM[iColumn]+lengthM[iColumn];j++) {
1452        index[nel]=indexM[j];
1453        element[nel++]=elementM[j];
1454      }
1455      which[n++]=iColumn;
1456      starts[n]=nel;
1457    }
1458    solution.deleteColumns(n,which);
1459    solution.dual();
1460    // Put back
1461    solution.addColumns(n,saveLower,saveUpper,saveObj,
1462                        starts,index,element);
1463    solution.dual();
1464    assert(eq(solution.objectiveValue(),1.5185098965e+03));
1465
1466    // reload with original
1467    solution.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
1468                         m.getObjCoefficients(),
1469                         m.getRowLower(),m.getRowUpper());
1470    // delete half rows
1471    n=0;
1472    nel=0;
1473    lower = m.getRowLower();
1474    upper = m.getRowUpper();
1475    startM = m.getMatrixByRow()->getVectorStarts();
1476    lengthM = m.getMatrixByRow()->getVectorLengths();
1477    indexM = m.getMatrixByRow()->getIndices();
1478    elementM = m.getMatrixByRow()->getElements();
1479    starts[0]=0;
1480    for (iRow=0;iRow<numberRows;iRow++) {
1481      if ((iRow&1)==0) {
1482        saveLower[n]=lower[iRow];
1483        saveUpper[n]=upper[iRow];
1484        int j;
1485        for (j=startM[iRow];j<startM[iRow]+lengthM[iRow];j++) {
1486          index[nel]=indexM[j];
1487          element[nel++]=elementM[j];
1488        }
1489        which[n++]=iRow;
1490        starts[n]=nel;
1491      }
1492    }
1493    solution.deleteRows(n,which);
1494    solution.dual();
1495    // Put back
1496    solution.addRows(n,saveLower,saveUpper,
1497                        starts,index,element);
1498    solution.dual();
1499    assert(eq(solution.objectiveValue(),1.5185098965e+03));
1500    solution.writeMps("yy.mps");
1501    // Delete all rows
1502    n=0;
1503    nel=0;
1504    lower = m.getRowLower();
1505    upper = m.getRowUpper();
1506    starts[0]=0;
1507    for (iRow=0;iRow<numberRows;iRow++) {
1508      saveLower[n]=lower[iRow];
1509      saveUpper[n]=upper[iRow];
1510      int j;
1511      for (j=startM[iRow];j<startM[iRow]+lengthM[iRow];j++) {
1512        index[nel]=indexM[j];
1513        element[nel++]=elementM[j];
1514      }
1515      which[n++]=iRow;
1516      starts[n]=nel;
1517    }
1518    solution.deleteRows(n,which);
1519    solution.dual();
1520    // Put back
1521    solution.addRows(n,saveLower,saveUpper,
1522                        starts,index,element);
1523    solution.dual();
1524    solution.writeMps("xx.mps");
1525    assert(eq(solution.objectiveValue(),1.5185098965e+03));
1526    // Zero out status array to give some interest
1527    memset(solution.statusArray()+numberColumns,0,numberRows);
1528    solution.primal(1);
1529    assert(eq(solution.objectiveValue(),1.5185098965e+03));
1530    // Delete all columns and rows
1531    n=0;
1532    for (iColumn=0;iColumn<numberColumns;iColumn++) {
1533      which[n++]=iColumn;
1534      starts[n]=nel;
1535    }
1536    solution.deleteColumns(n,which);
1537    n=0;
1538    for (iRow=0;iRow<numberRows;iRow++) {
1539      which[n++]=iRow;
1540      starts[n]=nel;
1541    }
1542    solution.deleteRows(n,which);
1543
1544    delete [] saveObj;
1545    delete [] saveLower;
1546    delete [] saveUpper;
1547    delete [] which;
1548    delete [] starts;
1549    delete [] index;
1550    delete [] element;
1551  }
1552#if 1
1553  // Test barrier
1554  {
1555    CoinMpsIO m;
1556    std::string fn = dirSample+"exmip1";
1557    m.readMps(fn.c_str(),"mps");
1558    ClpInterior solution;
1559    solution.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
1560                         m.getObjCoefficients(),
1561                         m.getRowLower(),m.getRowUpper());
1562    solution.primalDual();
1563  }
1564#endif
1565  // test network
1566#define QUADRATIC
1567  if (1) {   
1568    std::string fn = dirSample+"input.130";
1569    int numberColumns;
1570    int numberRows;
1571   
1572    FILE * fp = fopen(fn.c_str(),"r");
1573    if (!fp) {
1574      // Try in Data/Sample
1575      fn = "Data/Sample/input.130";
1576      fp = fopen(fn.c_str(),"r");
1577    }
1578    if (!fp) {
1579      fprintf(stderr,"Unable to open file input.130 in dirSample or Data/Sample directory\n");
1580    } else {
1581      int problem;
1582      char temp[100];
1583      // read and skip
1584      int x=fscanf(fp,"%s",temp);
1585      if (x<0)
1586        throw("bad fscanf");
1587      assert (!strcmp(temp,"BEGIN"));
1588      x=fscanf(fp,"%*s %*s %d %d %*s %*s %d %*s",&problem, &numberRows, 
1589             &numberColumns);
1590      if (x<0)
1591        throw("bad fscanf");
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 = static_cast<char>('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 = static_cast<char>('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.