source: branches/BSP/trunk/Clp/src/unitTest.cpp @ 1074

Last change on this file since 1074 was 1074, checked in by stefan, 12 years ago

error handling in case that an mps file from the netlib cannot be opened

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