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

Last change on this file since 1073 was 1073, checked in by ladanyi, 13 years ago

Got rid of dependency on Data/Netlib?

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 67.3 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        mps.readMps(fn.c_str(),"mps");
440        solutionBase.loadProblem(*mps.getMatrixByCol(),mps.getColLower(),
441                                 mps.getColUpper(),
442                                 mps.getObjCoefficients(),
443                                 mps.getRowLower(),mps.getRowUpper());
444       
445        solutionBase.setDblParam(ClpObjOffset,mps.objectiveOffset());
446      } 
447     
448      // Runs through strategies
449      if (algorithm==6||algorithm==7) {
450        // algorithms tested are at top of file
451        double testTime[NUMBER_ALGORITHMS];
452        std::string alg[NUMBER_ALGORITHMS];
453        int iTest;
454        for (iTest=0;iTest<NUMBER_ALGORITHMS;iTest++) {
455          ClpSolve solveOptions=setupForSolve(iTest,alg[iTest],1);
456          if (solveOptions.getSolveType()!=ClpSolve::notImplemented) {
457            double time1 = CoinCpuTime();
458            ClpSimplex solution=solutionBase;
459            if (solution.maximumSeconds()<0.0)
460              solution.setMaximumSeconds(120.0);
461            if (doVector) {
462              ClpMatrixBase * matrix = solution.clpMatrix();
463              if (dynamic_cast< ClpPackedMatrix*>(matrix)) {
464                ClpPackedMatrix * clpMatrix = dynamic_cast< ClpPackedMatrix*>(matrix);
465                clpMatrix->makeSpecialColumnCopy();
466              }
467            }
468            solution.initialSolve(solveOptions);
469            double time2 = CoinCpuTime()-time1;
470            testTime[iTest]=time2;
471            printf("Took %g seconds - status %d\n",time2,solution.problemStatus());
472            if (solution.problemStatus()) 
473              testTime[iTest]=1.0e20;
474          } else {
475            testTime[iTest]=1.0e30;
476          }
477        }
478        int iBest=-1;
479        double dBest=1.0e10;
480        printf("%s",fn.c_str());
481        for (iTest=0;iTest<NUMBER_ALGORITHMS;iTest++) {
482          if (testTime[iTest]<1.0e30) {
483            printf(" %s %g",
484                   alg[iTest].c_str(),testTime[iTest]);
485            if (testTime[iTest]<dBest) {
486              dBest=testTime[iTest];
487              iBest=iTest;
488            }
489          }
490        }
491        printf("\n");
492        if (iBest>=0)
493          printf("Best strategy for %s is %s (%d) which takes %g seconds\n",
494                 fn.c_str(),alg[iBest].c_str(),iBest,testTime[iBest]);
495        else
496          printf("No strategy finished in time limit\n");
497        continue;
498      }
499      double time1 = CoinCpuTime();
500      ClpSimplex solution=solutionBase;
501#if 0
502      solution.setOptimizationDirection(-1);
503      {
504        int j;
505        double * obj = solution.objective();
506        int n=solution.numberColumns();
507        for (j=0;j<n;j++)
508          obj[j] *= -1.0;
509      }
510#endif
511      ClpSolve::SolveType method;
512      ClpSolve::PresolveType presolveType;
513      ClpSolve solveOptions;
514      if (doPresolve)
515        presolveType=ClpSolve::presolveOn;
516      else
517        presolveType=ClpSolve::presolveOff;
518      solveOptions.setPresolveType(presolveType,5);
519      std::string nameAlgorithm;
520      if (algorithm!=5) {
521        if (algorithm==0) {
522          method=ClpSolve::useDual;
523          nameAlgorithm="dual";
524        } else if (algorithm==1) {
525          method=ClpSolve::usePrimalorSprint;
526          nameAlgorithm="primal";
527        } else if (algorithm==3) {
528          method=ClpSolve::automatic;
529          nameAlgorithm="either";
530        } else {
531          nameAlgorithm="barrier-slow";
532#ifdef UFL_BARRIER
533          solveOptions.setSpecialOption(4,4);
534          nameAlgorithm="barrier-UFL";
535#endif
536#ifdef WSSMP_BARRIER
537          solveOptions.setSpecialOption(4,2);
538          nameAlgorithm="barrier-WSSMP";
539#endif
540          method = ClpSolve::useBarrier;
541        }
542        solveOptions.setSolveType(method);
543      } else {
544        int iAlg = bestStrategy[m];
545        int presolveOff=iAlg/100;
546        iAlg=iAlg % 100;
547        if( !barrierAvailable&&iAlg==0) {
548          if (nRows[m]!=2172)
549            iAlg = 5; // try primal
550          else
551            iAlg=3; // d2q06c
552        }
553        solveOptions=setupForSolve(iAlg,nameAlgorithm,0);
554        if (presolveOff)
555          solveOptions.setPresolveType(ClpSolve::presolveOff);
556      }
557      if (doVector) {
558        ClpMatrixBase * matrix = solution.clpMatrix();
559        if (dynamic_cast< ClpPackedMatrix*>(matrix)) {
560          ClpPackedMatrix * clpMatrix = dynamic_cast< ClpPackedMatrix*>(matrix);
561          clpMatrix->makeSpecialColumnCopy();
562        }
563      }
564      solution.initialSolve(solveOptions);
565      double time2 = CoinCpuTime()-time1;
566      timeTaken += time2;
567      printf("%s took %g seconds using algorithm %s\n",fn.c_str(),time2,nameAlgorithm.c_str());
568      // Test objective solution value
569      {
570        double soln = solution.objectiveValue();
571        CoinRelFltEq eq(objValueTol[m]);
572        std::cerr <<soln <<",  " <<objValue[m] <<" diff "<<
573          soln-objValue[m]<<std::endl;
574        if(!eq(soln,objValue[m]))
575          printf("** difference fails\n");
576      }
577    }
578    printf("Total time %g seconds\n",timeTaken);
579  }
580 
581  testingMessage( "All tests completed successfully\n" );
582  return 0;
583}
584
585 
586// Display message on stdout and stderr
587void testingMessage( const char * const msg )
588{
589  std::cerr <<msg;
590  //cout <<endl <<"*****************************************"
591  //     <<endl <<msg <<endl;
592}
593
594//--------------------------------------------------------------------------
595// test factorization methods and simplex method and simple barrier
596void
597ClpSimplexUnitTest(const std::string & directory)
598{
599 
600  CoinRelFltEq eq(0.000001);
601
602  {
603    ClpSimplex solution;
604 
605    // matrix data
606    //deliberate hiccup of 2 between 0 and 1
607    CoinBigIndex start[5]={0,4,7,8,9};
608    int length[5]={2,3,1,1,1};
609    int rows[11]={0,2,-1,-1,0,1,2,0,1,2};
610    double elements[11]={7.0,2.0,1.0e10,1.0e10,-2.0,1.0,-2.0,1,1,1};
611    CoinPackedMatrix matrix(true,3,5,8,elements,rows,start,length);
612   
613    // rim data
614    double objective[7]={-4.0,1.0,0.0,0.0,0.0,0.0,0.0};
615    double rowLower[5]={14.0,3.0,3.0,1.0e10,1.0e10};
616    double rowUpper[5]={14.0,3.0,3.0,-1.0e10,-1.0e10};
617    double colLower[7]={0.0,0.0,0.0,0.0,0.0,0.0,0.0};
618    double colUpper[7]={100.0,100.0,100.0,100.0,100.0,100.0,100.0};
619   
620    // basis 1
621    int rowBasis1[3]={-1,-1,-1};
622    int colBasis1[5]={1,1,-1,-1,1};
623    solution.loadProblem(matrix,colLower,colUpper,objective,
624                         rowLower,rowUpper);
625    int i;
626    solution.createStatus();
627    for (i=0;i<3;i++) {
628      if (rowBasis1[i]<0) {
629        solution.setRowStatus(i,ClpSimplex::atLowerBound);
630      } else {
631        solution.setRowStatus(i,ClpSimplex::basic);
632      }
633    }
634    for (i=0;i<5;i++) {
635      if (colBasis1[i]<0) {
636        solution.setColumnStatus(i,ClpSimplex::atLowerBound);
637      } else {
638        solution.setColumnStatus(i,ClpSimplex::basic);
639      }
640    }
641    solution.setLogLevel(3+4+8+16+32);
642    solution.primal();
643    for (i=0;i<3;i++) {
644      if (rowBasis1[i]<0) {
645        solution.setRowStatus(i,ClpSimplex::atLowerBound);
646      } else {
647        solution.setRowStatus(i,ClpSimplex::basic);
648      }
649    }
650    for (i=0;i<5;i++) {
651      if (colBasis1[i]<0) {
652        solution.setColumnStatus(i,ClpSimplex::atLowerBound);
653      } else {
654        solution.setColumnStatus(i,ClpSimplex::basic);
655      }
656    }
657    // intricate stuff does not work with scaling
658    solution.scaling(0);
659    int returnCode = solution.factorize ( );
660    assert(!returnCode);
661    const double * colsol = solution.primalColumnSolution();
662    const double * rowsol = solution.primalRowSolution();
663    solution.getSolution(rowsol,colsol);
664    double colsol1[5]={20.0/7.0,3.0,0.0,0.0,23.0/7.0};
665    for (i=0;i<5;i++) {
666      assert(eq(colsol[i],colsol1[i]));
667    }
668    // now feed in again without actually doing factorization
669    ClpFactorization factorization2 = *solution.factorization();
670    ClpSimplex solution2 = solution;
671    solution2.setFactorization(factorization2);
672    solution2.createStatus();
673    for (i=0;i<3;i++) {
674      if (rowBasis1[i]<0) {
675        solution2.setRowStatus(i,ClpSimplex::atLowerBound);
676      } else {
677        solution2.setRowStatus(i,ClpSimplex::basic);
678      }
679    }
680    for (i=0;i<5;i++) {
681      if (colBasis1[i]<0) {
682        solution2.setColumnStatus(i,ClpSimplex::atLowerBound);
683      } else {
684        solution2.setColumnStatus(i,ClpSimplex::basic);
685      }
686    }
687    // intricate stuff does not work with scaling
688    solution2.scaling(0);
689    solution2.getSolution(rowsol,colsol);
690    colsol = solution2.primalColumnSolution();
691    rowsol = solution2.primalRowSolution();
692    for (i=0;i<5;i++) {
693      assert(eq(colsol[i],colsol1[i]));
694    }
695    solution2.setDualBound(0.1);
696    solution2.dual();
697    objective[2]=-1.0;
698    objective[3]=-0.5;
699    objective[4]=10.0;
700    solution.dual();
701    for (i=0;i<3;i++) {
702      rowLower[i]=-1.0e20;
703      colUpper[i+2]=0.0;
704    }
705    solution.setLogLevel(3);
706    solution.dual();
707    double rowObjective[]={1.0,0.5,-10.0};
708    solution.loadProblem(matrix,colLower,colUpper,objective,
709                         rowLower,rowUpper,rowObjective);
710    solution.dual();
711    solution.loadProblem(matrix,colLower,colUpper,objective,
712                         rowLower,rowUpper,rowObjective);
713    solution.primal();
714  }
715#ifndef COIN_NO_CLP_MESSAGE
716  {   
717    CoinMpsIO m;
718    std::string fn = directory+"exmip1";
719    m.readMps(fn.c_str(),"mps");
720    ClpSimplex solution;
721    solution.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
722                         m.getObjCoefficients(),
723                         m.getRowLower(),m.getRowUpper());
724    solution.dual();
725    // Test event handling
726    MyEventHandler handler;
727    solution.passInEventHandler(&handler);
728    int numberRows=solution.numberRows();
729    // make sure values pass has something to do
730    for (int i=0;i<numberRows;i++)
731      solution.setRowStatus(i,ClpSimplex::basic);
732    solution.primal(1);
733    assert (solution.secondaryStatus()==102); // Came out at end of pass
734  }
735  // Test Message handler
736  {   
737    CoinMpsIO m;
738    std::string fn = directory+"exmip1";
739    //fn = "Test/subGams4";
740    m.readMps(fn.c_str(),"mps");
741    ClpSimplex model;
742    model.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
743                         m.getObjCoefficients(),
744                         m.getRowLower(),m.getRowUpper());
745    // Message handler
746    MyMessageHandler messageHandler(&model);
747    std::cout<<"Testing derived message handler"<<std::endl;
748    model.passInMessageHandler(&messageHandler);
749    model.messagesPointer()->setDetailMessage(1,102);
750    model.setFactorizationFrequency(10);
751    model.primal();
752    model.primal(0,3);
753    model.setObjCoeff(3,-2.9473684210526314);
754    model.primal(0,3);
755    // Write saved solutions
756    int nc = model.getNumCols();
757    int s; 
758    std::deque<StdVectorDouble> fep = messageHandler.getFeasibleExtremePoints();
759    int numSavedSolutions = fep.size();
760    for ( s=0; s<numSavedSolutions; ++s ) {
761      const StdVectorDouble & solnVec = fep[s];
762      for ( int c=0; c<nc; ++c ) {
763        if (fabs(solnVec[c])>1.0e-8)
764          std::cout <<"Saved Solution: " <<s <<" ColNum: " <<c <<" Value: " <<solnVec[c] <<std::endl;
765      }
766    }
767    // Solve again without scaling
768    // and maximize then minimize
769    messageHandler.clearFeasibleExtremePoints();
770    model.scaling(0);
771    model.setOptimizationDirection(-1);
772    model.primal();
773    model.setOptimizationDirection(1);
774    model.primal();
775    fep = messageHandler.getFeasibleExtremePoints();
776    numSavedSolutions = fep.size();
777    for ( s=0; s<numSavedSolutions; ++s ) {
778      const StdVectorDouble & solnVec = fep[s];
779      for ( int c=0; c<nc; ++c ) {
780        if (fabs(solnVec[c])>1.0e-8)
781          std::cout <<"Saved Solution: " <<s <<" ColNum: " <<c <<" Value: " <<solnVec[c] <<std::endl;
782      }
783    }
784  }
785#endif
786  // Test dual ranging
787  {   
788    CoinMpsIO m;
789    std::string fn = directory+"exmip1";
790    m.readMps(fn.c_str(),"mps");
791    ClpSimplex model;
792    model.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
793                         m.getObjCoefficients(),
794                         m.getRowLower(),m.getRowUpper());
795    model.primal();
796    int which[13] = {0,1,2,3,4,5,6,7,8,9,10,11,12};
797    double costIncrease[13];
798    int sequenceIncrease[13];
799    double costDecrease[13];
800    int sequenceDecrease[13];
801    // ranging
802    model.dualRanging(13,which,costIncrease,sequenceIncrease,
803                      costDecrease,sequenceDecrease);
804    int i;
805    for ( i=0;i<13;i++)
806      printf("%d increase %g %d, decrease %g %d\n",
807             i,costIncrease[i],sequenceIncrease[i],
808             costDecrease[i],sequenceDecrease[i]);
809    assert (fabs(costDecrease[3])<1.0e-4);
810    assert (fabs(costIncrease[7]-1.0)<1.0e-4);
811    model.setOptimizationDirection(-1);
812    {
813      int j;
814      double * obj = model.objective();
815      int n=model.numberColumns();
816      for (j=0;j<n;j++) 
817        obj[j] *= -1.0;
818    }
819    double costIncrease2[13];
820    int sequenceIncrease2[13];
821    double costDecrease2[13];
822    int sequenceDecrease2[13];
823    // ranging
824    model.dualRanging(13,which,costIncrease2,sequenceIncrease2,
825                      costDecrease2,sequenceDecrease2);
826    for (i=0;i<13;i++) {
827      assert (fabs(costIncrease[i]-costDecrease2[i])<1.0e-6);
828      assert (fabs(costDecrease[i]-costIncrease2[i])<1.0e-6);
829      assert (sequenceIncrease[i]==sequenceDecrease2[i]);
830      assert (sequenceDecrease[i]==sequenceIncrease2[i]);
831    }
832    // Now delete all rows and see what happens
833    model.deleteRows(model.numberRows(),which);
834    model.primal();
835    // ranging
836    if (!model.dualRanging(8,which,costIncrease,sequenceIncrease,
837                           costDecrease,sequenceDecrease)) {
838      for (i=0;i<8;i++) {
839        printf("%d increase %g %d, decrease %g %d\n",
840               i,costIncrease[i],sequenceIncrease[i],
841               costDecrease[i],sequenceDecrease[i]);
842      }
843    }
844  }
845  // Test primal ranging
846  {   
847    CoinMpsIO m;
848    std::string fn = directory+"exmip1";
849    m.readMps(fn.c_str(),"mps");
850    ClpSimplex model;
851    model.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
852                         m.getObjCoefficients(),
853                         m.getRowLower(),m.getRowUpper());
854    model.primal();
855    int which[13] = {0,1,2,3,4,5,6,7,8,9,10,11,12};
856    double valueIncrease[13];
857    int sequenceIncrease[13];
858    double valueDecrease[13];
859    int sequenceDecrease[13];
860    // ranging
861    model.primalRanging(13,which,valueIncrease,sequenceIncrease,
862                      valueDecrease,sequenceDecrease);
863    int i;
864    for ( i=0;i<13;i++)
865      printf("%d increase %g %d, decrease %g %d\n",
866             i,valueIncrease[i],sequenceIncrease[i],
867             valueDecrease[i],sequenceDecrease[i]);
868    assert (fabs(valueIncrease[3]-0.642857)<1.0e-4);
869    assert (fabs(valueIncrease[8]-2.95113)<1.0e-4);
870#if 0
871    // out until I find optimization bug
872    // Test parametrics
873    ClpSimplexOther * model2 = (ClpSimplexOther *) (&model);
874    double rhs[]={ 1.0,2.0,3.0,4.0,5.0};
875    double endingTheta=1.0;
876    model2->scaling(0);
877    model2->setLogLevel(63);
878    model2->parametrics(0.0,endingTheta,0.1,
879                        NULL,NULL,rhs,rhs,NULL);
880#endif
881  }
882  // Test binv etc
883  {   
884    /*
885       Wolsey : Page 130
886       max 4x1 -  x2
887       7x1 - 2x2    <= 14
888       x2    <= 3
889       2x1 - 2x2    <= 3
890       x1 in Z+, x2 >= 0
891
892       note slacks are -1 in Clp so signs may be different
893    */
894   
895    int n_cols = 2;
896    int n_rows = 3;
897   
898    double obj[2] = {-4.0, 1.0};
899    double collb[2] = {0.0, 0.0};
900    double colub[2] = {COIN_DBL_MAX, COIN_DBL_MAX};
901    double rowlb[3] = {-COIN_DBL_MAX, -COIN_DBL_MAX, -COIN_DBL_MAX};
902    double rowub[3] = {14.0, 3.0, 3.0};
903   
904    int rowIndices[5] =  {0,     2,    0,    1,    2};
905    int colIndices[5] =  {0,     0,    1,    1,    1};
906    double elements[5] = {7.0, 2.0, -2.0,  1.0, -2.0};
907    CoinPackedMatrix M(true, rowIndices, colIndices, elements, 5);
908
909    ClpSimplex model;
910    model.loadProblem(M, collb, colub, obj, rowlb, rowub);
911    model.dual(0,1); // keep factorization
912   
913    //check that the tableau matches wolsey (B-1 A)
914    // slacks in second part of binvA
915    double * binvA = (double*) malloc((n_cols+n_rows) * sizeof(double));
916   
917    printf("B-1 A by row\n");
918    int i;
919    for( i = 0; i < n_rows; i++){
920      model.getBInvARow(i, binvA,binvA+n_cols);
921      printf("row: %d -> ",i);
922      for(int j=0; j < n_cols+n_rows; j++){
923        printf("%g, ", binvA[j]);
924      }
925      printf("\n");
926    }
927    // See if can re-use factorization AND arrays
928    model.primal(0,3+4); // keep factorization
929    // And do by column
930    printf("B-1 A by column\n");
931    for( i = 0; i < n_rows+n_cols; i++){
932      model.getBInvACol(i, binvA);
933      printf("column: %d -> ",i);
934      for(int j=0; j < n_rows; j++){
935        printf("%g, ", binvA[j]);
936      }
937      printf("\n");
938    }
939    free(binvA);
940    model.setColUpper(1,2.0);
941    model.dual(0,2+4); // use factorization and arrays
942    model.dual(0,2); // hopefully will not use factorization
943    model.primal(0,3+4); // keep factorization
944    // but say basis has changed
945    model.setWhatsChanged(model.whatsChanged()&(~512));
946    model.dual(0,2); // hopefully will not use factorization
947  }
948  // test steepest edge
949  {   
950    CoinMpsIO m;
951    std::string fn = directory+"finnis";
952    int returnCode = m.readMps(fn.c_str(),"mps");
953    if (returnCode) {
954      // probable cause is that gz not there
955      fprintf(stderr,"Unable to open finnis.mps in %s!\n", directory.c_str());
956      fprintf(stderr,"Most probable cause is finnis.mps is gzipped i.e. finnis.mps.gz and libz has not been activated\n");
957      fprintf(stderr,"Either gunzip files or edit Makefiles/Makefile.location to get libz\n");
958      exit(999);
959    }
960    ClpModel model;
961    model.loadProblem(*m.getMatrixByCol(),m.getColLower(),
962                    m.getColUpper(),
963                    m.getObjCoefficients(),
964                    m.getRowLower(),m.getRowUpper());
965    ClpSimplex solution(model);
966
967    solution.scaling(1); 
968    solution.setDualBound(1.0e8);
969    //solution.factorization()->maximumPivots(1);
970    //solution.setLogLevel(3);
971    solution.setDualTolerance(1.0e-7);
972    // set objective sense,
973    ClpDualRowSteepest steep;
974    solution.setDualRowPivotAlgorithm(steep);
975    solution.setDblParam(ClpObjOffset,m.objectiveOffset());
976    solution.dual();
977  }
978  // test normal solution
979  {   
980    CoinMpsIO m;
981    std::string fn = directory+"afiro";
982    m.readMps(fn.c_str(),"mps");
983    ClpSimplex solution;
984    ClpModel model;
985    // do twice - without and with scaling
986    int iPass;
987    for (iPass=0;iPass<2;iPass++) {
988      // explicit row objective for testing
989      int nr = m.getNumRows();
990      double * rowObj = new double[nr];
991      CoinFillN(rowObj,nr,0.0);
992      model.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
993                      m.getObjCoefficients(),
994                      m.getRowLower(),m.getRowUpper(),rowObj);
995      delete [] rowObj;
996      solution = ClpSimplex(model);
997      if (iPass) {
998        solution.scaling();
999      }
1000      solution.dual();
1001      solution.dual();
1002      // test optimal
1003      assert (solution.status()==0);
1004      int numberColumns = solution.numberColumns();
1005      int numberRows = solution.numberRows();
1006      CoinPackedVector colsol(numberColumns,solution.primalColumnSolution());
1007      double * objective = solution.objective();
1008      double objValue = colsol.dotProduct(objective);
1009      CoinRelFltEq eq(1.0e-8);
1010      assert(eq(objValue,-4.6475314286e+02));
1011      // Test auxiliary model
1012      //solution.scaling(0);
1013      solution.auxiliaryModel(63-2); // bounds may change
1014      solution.dual();
1015      solution.primal();
1016      solution.allSlackBasis();
1017      solution.dual();
1018      assert(eq(solution.objectiveValue(),-4.6475314286e+02));
1019      solution.auxiliaryModel(-1);
1020      solution.dual();
1021      assert(eq(solution.objectiveValue(),-4.6475314286e+02));
1022      double * lower = solution.columnLower();
1023      double * upper = solution.columnUpper();
1024      double * sol = solution.primalColumnSolution();
1025      double * result = new double[numberColumns];
1026      CoinFillN ( result, numberColumns,0.0);
1027      solution.matrix()->transposeTimes(solution.dualRowSolution(), result);
1028      int iRow , iColumn;
1029      // see if feasible and dual feasible
1030      for (iColumn=0;iColumn<numberColumns;iColumn++) {
1031        double value = sol[iColumn];
1032        assert(value<upper[iColumn]+1.0e-8);
1033        assert(value>lower[iColumn]-1.0e-8);
1034        value = objective[iColumn]-result[iColumn];
1035        assert (value>-1.0e-5);
1036        if (sol[iColumn]>1.0e-5)
1037          assert (value<1.0e-5);
1038      }
1039      delete [] result;
1040      result = new double[numberRows];
1041      CoinFillN ( result, numberRows,0.0);
1042      solution.matrix()->times(colsol, result);
1043      lower = solution.rowLower();
1044      upper = solution.rowUpper();
1045      sol = solution.primalRowSolution();
1046      for (iRow=0;iRow<numberRows;iRow++) {
1047        double value = result[iRow];
1048        assert(eq(value,sol[iRow]));
1049        assert(value<upper[iRow]+1.0e-8);
1050        assert(value>lower[iRow]-1.0e-8);
1051      }
1052      delete [] result;
1053      // test row objective
1054      double * rowObjective = solution.rowObjective();
1055      CoinDisjointCopyN(solution.dualRowSolution(),numberRows,rowObjective);
1056      CoinDisjointCopyN(solution.dualColumnSolution(),numberColumns,objective);
1057      // this sets up all slack basis
1058      solution.createStatus();
1059      solution.dual();
1060      CoinFillN(rowObjective,numberRows,0.0);
1061      CoinDisjointCopyN(m.getObjCoefficients(),numberColumns,objective);
1062      solution.dual();
1063    }
1064  }
1065  // test unbounded
1066  {   
1067    CoinMpsIO m;
1068    std::string fn = directory+"brandy";
1069    m.readMps(fn.c_str(),"mps");
1070    ClpSimplex solution;
1071    // do twice - without and with scaling
1072    int iPass;
1073    for (iPass=0;iPass<2;iPass++) {
1074      solution.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
1075                      m.getObjCoefficients(),
1076                      m.getRowLower(),m.getRowUpper());
1077      if (iPass)
1078        solution.scaling();
1079      solution.setOptimizationDirection(-1);
1080      // test unbounded and ray
1081#ifdef DUAL
1082      solution.setDualBound(100.0);
1083      solution.dual();
1084#else
1085      solution.primal();
1086#endif
1087      assert (solution.status()==2);
1088      int numberColumns = solution.numberColumns();
1089      int numberRows = solution.numberRows();
1090      double * lower = solution.columnLower();
1091      double * upper = solution.columnUpper();
1092      double * sol = solution.primalColumnSolution();
1093      double * ray = solution.unboundedRay();
1094      double * objective = solution.objective();
1095      double objChange=0.0;
1096      int iRow , iColumn;
1097      // make sure feasible and columns form ray
1098      for (iColumn=0;iColumn<numberColumns;iColumn++) {
1099        double value = sol[iColumn];
1100        assert(value<upper[iColumn]+1.0e-8);
1101        assert(value>lower[iColumn]-1.0e-8);
1102        value = ray[iColumn];
1103        if (value>0.0)
1104          assert(upper[iColumn]>1.0e30);
1105        else if (value<0.0)
1106          assert(lower[iColumn]<-1.0e30);
1107        objChange += value*objective[iColumn];
1108      }
1109      // make sure increasing objective
1110      assert(objChange>0.0);
1111      double * result = new double[numberRows];
1112      CoinFillN ( result, numberRows,0.0);
1113      solution.matrix()->times(sol, result);
1114      lower = solution.rowLower();
1115      upper = solution.rowUpper();
1116      sol = solution.primalRowSolution();
1117      for (iRow=0;iRow<numberRows;iRow++) {
1118        double value = result[iRow];
1119        assert(eq(value,sol[iRow]));
1120        assert(value<upper[iRow]+2.0e-8);
1121        assert(value>lower[iRow]-2.0e-8);
1122      }
1123      CoinFillN ( result, numberRows,0.0);
1124      solution.matrix()->times(ray, result);
1125      // there may be small differences (especially if scaled)
1126      for (iRow=0;iRow<numberRows;iRow++) {
1127        double value = result[iRow];
1128        if (value>1.0e-8)
1129          assert(upper[iRow]>1.0e30);
1130        else if (value<-1.0e-8)
1131          assert(lower[iRow]<-1.0e30);
1132      }
1133      delete [] result;
1134      delete [] ray;
1135    }
1136  }
1137  // test infeasible
1138  {   
1139    CoinMpsIO m;
1140    std::string fn = directory+"brandy";
1141    m.readMps(fn.c_str(),"mps");
1142    ClpSimplex solution;
1143    // do twice - without and with scaling
1144    int iPass;
1145    for (iPass=0;iPass<2;iPass++) {
1146      solution.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
1147                      m.getObjCoefficients(),
1148                      m.getRowLower(),m.getRowUpper());
1149      if (iPass)
1150        solution.scaling();
1151      // test infeasible and ray
1152      solution.columnUpper()[0]=0.0;
1153#ifdef DUAL
1154      solution.setDualBound(100.0);
1155      solution.dual();
1156#else
1157      solution.primal();
1158#endif
1159      assert (solution.status()==1);
1160      int numberColumns = solution.numberColumns();
1161      int numberRows = solution.numberRows();
1162      double * lower = solution.rowLower();
1163      double * upper = solution.rowUpper();
1164      double * ray = solution.infeasibilityRay();
1165      assert(ray);
1166      // construct proof of infeasibility
1167      int iRow , iColumn;
1168      double lo=0.0,up=0.0;
1169      int nl=0,nu=0;
1170      for (iRow=0;iRow<numberRows;iRow++) {
1171        if (lower[iRow]>-1.0e20) {
1172          lo += ray[iRow]*lower[iRow];
1173        } else {
1174          if (ray[iRow]>1.0e-8) 
1175            nl++;
1176        }
1177        if (upper[iRow]<1.0e20) {
1178          up += ray[iRow]*upper[iRow];
1179        } else {
1180          if (ray[iRow]>1.0e-8) 
1181            nu++;
1182        }
1183      }
1184      if (nl)
1185        lo=-1.0e100;
1186      if (nu)
1187        up=1.0e100;
1188      double * result = new double[numberColumns];
1189      double lo2=0.0,up2=0.0;
1190      CoinFillN ( result, numberColumns,0.0);
1191      solution.matrix()->transposeTimes(ray, result);
1192      lower = solution.columnLower();
1193      upper = solution.columnUpper();
1194      nl=nu=0;
1195      for (iColumn=0;iColumn<numberColumns;iColumn++) {
1196        if (result[iColumn]>1.0e-8) {
1197          if (lower[iColumn]>-1.0e20)
1198            lo2 += result[iColumn]*lower[iColumn];
1199          else
1200            nl++;
1201          if (upper[iColumn]<1.0e20)
1202            up2 += result[iColumn]*upper[iColumn];
1203          else
1204            nu++;
1205        } else if (result[iColumn]<-1.0e-8) {
1206          if (lower[iColumn]>-1.0e20)
1207            up2 += result[iColumn]*lower[iColumn];
1208          else
1209            nu++;
1210          if (upper[iColumn]<1.0e20)
1211            lo2 += result[iColumn]*upper[iColumn];
1212          else
1213            nl++;
1214        }
1215      }
1216      if (nl)
1217        lo2=-1.0e100;
1218      if (nu)
1219        up2=1.0e100;
1220      // make sure inconsistency
1221      assert(lo2>up||up2<lo);
1222      delete [] result;
1223      delete [] ray;
1224    }
1225  }
1226  // test delete and add
1227  {   
1228    CoinMpsIO m;
1229    std::string fn = directory+"brandy";
1230    m.readMps(fn.c_str(),"mps");
1231    ClpSimplex solution;
1232    solution.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
1233                         m.getObjCoefficients(),
1234                         m.getRowLower(),m.getRowUpper());
1235    solution.dual();
1236    CoinRelFltEq eq(1.0e-8);
1237    assert(eq(solution.objectiveValue(),1.5185098965e+03));
1238
1239    int numberColumns = solution.numberColumns();
1240    int numberRows = solution.numberRows();
1241    double * saveObj = new double [numberColumns];
1242    double * saveLower = new double[numberRows+numberColumns];
1243    double * saveUpper = new double[numberRows+numberColumns];
1244    int * which = new int [numberRows+numberColumns];
1245
1246    int numberElements = m.getMatrixByCol()->getNumElements();
1247    int * starts = new int[numberRows+numberColumns];
1248    int * index = new int[numberElements];
1249    double * element = new double[numberElements];
1250
1251    const CoinBigIndex * startM;
1252    const int * lengthM;
1253    const int * indexM;
1254    const double * elementM;
1255
1256    int n,nel;
1257
1258    // delete non basic columns
1259    n=0;
1260    nel=0;
1261    int iRow , iColumn;
1262    const double * lower = m.getColLower();
1263    const double * upper = m.getColUpper();
1264    const double * objective = m.getObjCoefficients();
1265    startM = m.getMatrixByCol()->getVectorStarts();
1266    lengthM = m.getMatrixByCol()->getVectorLengths();
1267    indexM = m.getMatrixByCol()->getIndices();
1268    elementM = m.getMatrixByCol()->getElements();
1269    starts[0]=0;
1270    for (iColumn=0;iColumn<numberColumns;iColumn++) {
1271      if (solution.getColumnStatus(iColumn)!=ClpSimplex::basic) {
1272        saveObj[n]=objective[iColumn];
1273        saveLower[n]=lower[iColumn];
1274        saveUpper[n]=upper[iColumn];
1275        int j;
1276        for (j=startM[iColumn];j<startM[iColumn]+lengthM[iColumn];j++) {
1277          index[nel]=indexM[j];
1278          element[nel++]=elementM[j];
1279        }
1280        which[n++]=iColumn;
1281        starts[n]=nel;
1282      }
1283    }
1284    solution.deleteColumns(n,which);
1285    solution.dual();
1286    // Put back
1287    solution.addColumns(n,saveLower,saveUpper,saveObj,
1288                        starts,index,element);
1289    solution.dual();
1290    assert(eq(solution.objectiveValue(),1.5185098965e+03));
1291    // Delete all columns and add back
1292    n=0;
1293    nel=0;
1294    starts[0]=0;
1295    lower = m.getColLower();
1296    upper = m.getColUpper();
1297    objective = m.getObjCoefficients();
1298    for (iColumn=0;iColumn<numberColumns;iColumn++) {
1299      saveObj[n]=objective[iColumn];
1300      saveLower[n]=lower[iColumn];
1301      saveUpper[n]=upper[iColumn];
1302      int j;
1303      for (j=startM[iColumn];j<startM[iColumn]+lengthM[iColumn];j++) {
1304        index[nel]=indexM[j];
1305        element[nel++]=elementM[j];
1306      }
1307      which[n++]=iColumn;
1308      starts[n]=nel;
1309    }
1310    solution.deleteColumns(n,which);
1311    solution.dual();
1312    // Put back
1313    solution.addColumns(n,saveLower,saveUpper,saveObj,
1314                        starts,index,element);
1315    solution.dual();
1316    assert(eq(solution.objectiveValue(),1.5185098965e+03));
1317
1318    // reload with original
1319    solution.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
1320                         m.getObjCoefficients(),
1321                         m.getRowLower(),m.getRowUpper());
1322    // delete half rows
1323    n=0;
1324    nel=0;
1325    lower = m.getRowLower();
1326    upper = m.getRowUpper();
1327    startM = m.getMatrixByRow()->getVectorStarts();
1328    lengthM = m.getMatrixByRow()->getVectorLengths();
1329    indexM = m.getMatrixByRow()->getIndices();
1330    elementM = m.getMatrixByRow()->getElements();
1331    starts[0]=0;
1332    for (iRow=0;iRow<numberRows;iRow++) {
1333      if ((iRow&1)==0) {
1334        saveLower[n]=lower[iRow];
1335        saveUpper[n]=upper[iRow];
1336        int j;
1337        for (j=startM[iRow];j<startM[iRow]+lengthM[iRow];j++) {
1338          index[nel]=indexM[j];
1339          element[nel++]=elementM[j];
1340        }
1341        which[n++]=iRow;
1342        starts[n]=nel;
1343      }
1344    }
1345    solution.deleteRows(n,which);
1346    solution.dual();
1347    // Put back
1348    solution.addRows(n,saveLower,saveUpper,
1349                        starts,index,element);
1350    solution.dual();
1351    assert(eq(solution.objectiveValue(),1.5185098965e+03));
1352    solution.writeMps("yy.mps");
1353    // Delete all rows
1354    n=0;
1355    nel=0;
1356    lower = m.getRowLower();
1357    upper = m.getRowUpper();
1358    starts[0]=0;
1359    for (iRow=0;iRow<numberRows;iRow++) {
1360      saveLower[n]=lower[iRow];
1361      saveUpper[n]=upper[iRow];
1362      int j;
1363      for (j=startM[iRow];j<startM[iRow]+lengthM[iRow];j++) {
1364        index[nel]=indexM[j];
1365        element[nel++]=elementM[j];
1366      }
1367      which[n++]=iRow;
1368      starts[n]=nel;
1369    }
1370    solution.deleteRows(n,which);
1371    solution.dual();
1372    // Put back
1373    solution.addRows(n,saveLower,saveUpper,
1374                        starts,index,element);
1375    solution.dual();
1376    solution.writeMps("xx.mps");
1377    assert(eq(solution.objectiveValue(),1.5185098965e+03));
1378    // Zero out status array to give some interest
1379    memset(solution.statusArray()+numberColumns,0,numberRows);
1380    solution.primal(1);
1381    assert(eq(solution.objectiveValue(),1.5185098965e+03));
1382    // Delete all columns and rows
1383    n=0;
1384    for (iColumn=0;iColumn<numberColumns;iColumn++) {
1385      which[n++]=iColumn;
1386      starts[n]=nel;
1387    }
1388    solution.deleteColumns(n,which);
1389    n=0;
1390    for (iRow=0;iRow<numberRows;iRow++) {
1391      which[n++]=iRow;
1392      starts[n]=nel;
1393    }
1394    solution.deleteRows(n,which);
1395
1396    delete [] saveObj;
1397    delete [] saveLower;
1398    delete [] saveUpper;
1399    delete [] which;
1400    delete [] starts;
1401    delete [] index;
1402    delete [] element;
1403  }
1404#if 1
1405  // Test barrier
1406  {
1407    CoinMpsIO m;
1408    std::string fn = directory+"exmip1";
1409    m.readMps(fn.c_str(),"mps");
1410    ClpInterior solution;
1411    solution.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
1412                         m.getObjCoefficients(),
1413                         m.getRowLower(),m.getRowUpper());
1414    solution.primalDual();
1415  }
1416#endif
1417  // test network
1418#define QUADRATIC
1419  if (1) {   
1420    std::string fn = directory+"input.130";
1421    int numberColumns;
1422    int numberRows;
1423   
1424    FILE * fp = fopen(fn.c_str(),"r");
1425    if (!fp) {
1426      fprintf(stderr,"Unable to open file input.130 in directory %s\n",
1427              directory.c_str());
1428    } else {
1429      int problem;
1430      char temp[100];
1431      // read and skip
1432      fscanf(fp,"%s",temp);
1433      assert (!strcmp(temp,"BEGIN"));
1434      fscanf(fp,"%*s %*s %d %d %*s %*s %d %*s",&problem, &numberRows, 
1435             &numberColumns);
1436      // scan down to SUPPLY
1437      while (fgets(temp,100,fp)) {
1438        if (!strncmp(temp,"SUPPLY",6))
1439          break;
1440      }
1441      if (strncmp(temp,"SUPPLY",6)) {
1442        fprintf(stderr,"Unable to find SUPPLY\n");
1443        exit(2);
1444      }
1445      // get space for rhs
1446      double * lower = new double[numberRows];
1447      double * upper = new double[numberRows];
1448      int i;
1449      for (i=0;i<numberRows;i++) {
1450        lower[i]=0.0;
1451        upper[i]=0.0;
1452      }
1453      // ***** Remember to convert to C notation
1454      while (fgets(temp,100,fp)) {
1455        int row;
1456        int value;
1457        if (!strncmp(temp,"ARCS",4))
1458          break;
1459        sscanf(temp,"%d %d",&row,&value);
1460        upper[row-1]=-value;
1461        lower[row-1]=-value;
1462      }
1463      if (strncmp(temp,"ARCS",4)) {
1464        fprintf(stderr,"Unable to find ARCS\n");
1465        exit(2);
1466      }
1467      // number of columns may be underestimate
1468      int * head = new int[2*numberColumns];
1469      int * tail = new int[2*numberColumns];
1470      int * ub = new int[2*numberColumns];
1471      int * cost = new int[2*numberColumns];
1472      // ***** Remember to convert to C notation
1473      numberColumns=0;
1474      while (fgets(temp,100,fp)) {
1475        int iHead;
1476        int iTail;
1477        int iUb;
1478        int iCost;
1479        if (!strncmp(temp,"DEMAND",6))
1480          break;
1481        sscanf(temp,"%d %d %d %d",&iHead,&iTail,&iCost,&iUb);
1482        iHead--;
1483        iTail--;
1484        head[numberColumns]=iHead;
1485        tail[numberColumns]=iTail;
1486        ub[numberColumns]=iUb;
1487        cost[numberColumns]=iCost;
1488        numberColumns++;
1489      }
1490      if (strncmp(temp,"DEMAND",6)) {
1491        fprintf(stderr,"Unable to find DEMAND\n");
1492        exit(2);
1493      }
1494      // ***** Remember to convert to C notation
1495      while (fgets(temp,100,fp)) {
1496        int row;
1497        int value;
1498        if (!strncmp(temp,"END",3))
1499          break;
1500        sscanf(temp,"%d %d",&row,&value);
1501        upper[row-1]=value;
1502        lower[row-1]=value;
1503      }
1504      if (strncmp(temp,"END",3)) {
1505        fprintf(stderr,"Unable to find END\n");
1506        exit(2);
1507      }
1508      printf("Problem %d has %d rows and %d columns\n",problem,
1509             numberRows,numberColumns);
1510      fclose(fp);
1511      ClpSimplex  model;
1512      // now build model
1513     
1514      double * objective =new double[numberColumns];
1515      double * lowerColumn = new double[numberColumns];
1516      double * upperColumn = new double[numberColumns];
1517     
1518      double * element = new double [2*numberColumns];
1519      CoinBigIndex * start = new CoinBigIndex [numberColumns+1];
1520      int * row = new int[2*numberColumns];
1521      start[numberColumns]=2*numberColumns;
1522      for (i=0;i<numberColumns;i++) {
1523        start[i]=2*i;
1524        element[2*i]=-1.0;
1525        element[2*i+1]=1.0;
1526        row[2*i]=head[i];
1527        row[2*i+1]=tail[i];
1528        lowerColumn[i]=0.0;
1529        upperColumn[i]=ub[i];
1530        objective[i]=cost[i];
1531      }
1532      // Create Packed Matrix
1533      CoinPackedMatrix matrix;
1534      int * lengths = NULL;
1535      matrix.assignMatrix(true,numberRows,numberColumns,
1536                          2*numberColumns,element,row,start,lengths);
1537      // load model
1538      model.loadProblem(matrix,
1539                        lowerColumn,upperColumn,objective,
1540                        lower,upper);
1541      model.factorization()->maximumPivots(200+model.numberRows()/100);
1542      model.createStatus();
1543      double time1 = CoinCpuTime();
1544      model.dual();
1545      std::cout<<"Network problem, ClpPackedMatrix took "<<CoinCpuTime()-time1<<" seconds"<<std::endl;
1546      ClpPlusMinusOneMatrix * plusMinus = new ClpPlusMinusOneMatrix(matrix);
1547      assert (plusMinus->getIndices()); // would be zero if not +- one
1548      //ClpPlusMinusOneMatrix *plusminus_matrix;
1549
1550      //plusminus_matrix = new ClpPlusMinusOneMatrix;
1551
1552      //plusminus_matrix->passInCopy(numberRows, numberColumns, true, plusMinus->getMutableIndices(),
1553      //                         plusMinus->startPositive(),plusMinus->startNegative());
1554      model.loadProblem(*plusMinus,
1555                        lowerColumn,upperColumn,objective,
1556                        lower,upper);
1557      //model.replaceMatrix( plusminus_matrix , true);
1558      delete plusMinus;
1559      //model.createStatus();
1560      //model.initialSolve();
1561      //model.writeMps("xx.mps");
1562     
1563      model.factorization()->maximumPivots(200+model.numberRows()/100);
1564      model.createStatus();
1565      time1 = CoinCpuTime();
1566      model.dual();
1567      std::cout<<"Network problem, ClpPlusMinusOneMatrix took "<<CoinCpuTime()-time1<<" seconds"<<std::endl;
1568      ClpNetworkMatrix network(numberColumns,head,tail);
1569      model.loadProblem(network,
1570                        lowerColumn,upperColumn,objective,
1571                        lower,upper);
1572     
1573      model.factorization()->maximumPivots(200+model.numberRows()/100);
1574      model.createStatus();
1575      time1 = CoinCpuTime();
1576      model.dual();
1577      std::cout<<"Network problem, ClpNetworkMatrix took "<<CoinCpuTime()-time1<<" seconds"<<std::endl;
1578      delete [] lower;
1579      delete [] upper;
1580      delete [] head;
1581      delete [] tail;
1582      delete [] ub;
1583      delete [] cost;
1584      delete [] objective;
1585      delete [] lowerColumn;
1586      delete [] upperColumn;
1587    }
1588  }
1589#ifdef QUADRATIC
1590  // Test quadratic to solve linear
1591  if (1) {   
1592    CoinMpsIO m;
1593    std::string fn = directory+"exmip1";
1594    m.readMps(fn.c_str(),"mps");
1595    ClpSimplex solution;
1596    solution.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
1597                         m.getObjCoefficients(),
1598                         m.getRowLower(),m.getRowUpper());
1599    //solution.dual();
1600    // get quadratic part
1601    int numberColumns=solution.numberColumns();
1602    int * start=new int [numberColumns+1];
1603    int * column = new int[numberColumns];
1604    double * element = new double[numberColumns];
1605    int i;
1606    start[0]=0;
1607    int n=0;
1608    int kk=numberColumns-1;
1609    int kk2=numberColumns-1;
1610    for (i=0;i<numberColumns;i++) {
1611      if (i>=kk) {
1612        column[n]=i;
1613        if (i>=kk2)
1614          element[n]=1.0e-1;
1615        else
1616          element[n]=0.0;
1617        n++;
1618      }
1619      start[i+1]=n;
1620    }
1621    // Load up objective
1622    solution.loadQuadraticObjective(numberColumns,start,column,element);
1623    delete [] start;
1624    delete [] column;
1625    delete [] element;
1626    //solution.quadraticSLP(50,1.0e-4);
1627    CoinRelFltEq eq(1.0e-4);
1628    //assert(eq(objValue,820.0));
1629    //solution.setLogLevel(63);
1630    solution.primal();
1631    printSol(solution);
1632    //assert(eq(objValue,3.2368421));
1633    //exit(77);
1634  }
1635  // Test quadratic
1636  if (1) {   
1637    CoinMpsIO m;
1638    std::string fn = directory+"share2qp";
1639    //fn = "share2qpb";
1640    m.readMps(fn.c_str(),"mps");
1641    ClpSimplex model;
1642    model.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
1643                         m.getObjCoefficients(),
1644                         m.getRowLower(),m.getRowUpper());
1645    model.dual();
1646    // get quadratic part
1647    int * start=NULL;
1648    int * column = NULL;
1649    double * element = NULL;
1650    m.readQuadraticMps(NULL,start,column,element,2);
1651    int column2[200];
1652    double element2[200];
1653    int start2[80];
1654    int j;
1655    start2[0]=0;
1656    int nel=0;
1657    bool good=false;
1658    for (j=0;j<79;j++) {
1659      if (start[j]==start[j+1]) {
1660        column2[nel]=j;
1661        element2[nel]=0.0;
1662        nel++;
1663      } else {
1664        int i;
1665        for (i=start[j];i<start[j+1];i++) {
1666          column2[nel]=column[i];
1667          element2[nel++]=element[i];
1668        }
1669      }
1670      start2[j+1]=nel;
1671    }
1672    // Load up objective
1673    if (good)
1674      model.loadQuadraticObjective(model.numberColumns(),start2,column2,element2);
1675    else
1676      model.loadQuadraticObjective(model.numberColumns(),start,column,element);
1677    delete [] start;
1678    delete [] column;
1679    delete [] element;
1680    int numberColumns=model.numberColumns();
1681    model.scaling(0);
1682#if 0
1683    model.nonlinearSLP(50,1.0e-4);
1684#else
1685    // Get feasible
1686    ClpObjective * saveObjective = model.objectiveAsObject()->clone();
1687    ClpLinearObjective zeroObjective(NULL,numberColumns);
1688    model.setObjective(&zeroObjective);
1689    model.dual();
1690    model.setObjective(saveObjective);
1691    delete saveObjective;
1692#endif
1693    //model.setLogLevel(63);
1694    //exit(77);
1695    model.setFactorizationFrequency(10);
1696    model.primal();
1697    printSol(model);
1698    double objValue = model.getObjValue();
1699    CoinRelFltEq eq(1.0e-4);
1700    assert(eq(objValue,-400.92));
1701    // and again for barrier
1702    model.barrier(false);
1703    //printSol(model);
1704    model.allSlackBasis();
1705    model.primal();
1706    //printSol(model);
1707  }
1708  if (0) {   
1709    CoinMpsIO m;
1710    std::string fn = "./beale";
1711    //fn = "./jensen";
1712    m.readMps(fn.c_str(),"mps");
1713    ClpSimplex solution;
1714    solution.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
1715                         m.getObjCoefficients(),
1716                         m.getRowLower(),m.getRowUpper());
1717    solution.setDblParam(ClpObjOffset,m.objectiveOffset());
1718    solution.dual();
1719    // get quadratic part
1720    int * start=NULL;
1721    int * column = NULL;
1722    double * element = NULL;
1723    m.readQuadraticMps(NULL,start,column,element,2);
1724    // Load up objective
1725    solution.loadQuadraticObjective(solution.numberColumns(),start,column,element);
1726    delete [] start;
1727    delete [] column;
1728    delete [] element;
1729    solution.primal(1);
1730    solution.nonlinearSLP(50,1.0e-4);
1731    double objValue = solution.getObjValue();
1732    CoinRelFltEq eq(1.0e-4);
1733    assert(eq(objValue,0.5));
1734    solution.primal();
1735    objValue = solution.getObjValue();
1736    assert(eq(objValue,0.5));
1737  }
1738#endif 
1739}
Note: See TracBrowser for help on using the repository browser.