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

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

get rid of compiler warnings

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