source: branches/pre/Test/unitTest.cpp @ 186

Last change on this file since 186 was 186, checked in by forrest, 18 years ago

stuff

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 47.9 KB
Line 
1// Copyright (C) 2002, International Business Machines
2// Corporation and others.  All Rights Reserved.
3
4#include "CoinPragma.hpp"
5#include <cassert>
6#include <cstdio>
7#include <cmath>
8#include <cfloat>
9#include <string>
10#include <iostream>
11
12#include "CoinMpsIO.hpp"
13#include "CoinPackedMatrix.hpp"
14#include "CoinPackedVector.hpp"
15#include "CoinHelperFunctions.hpp"
16
17#include "ClpFactorization.hpp"
18#include "ClpSimplex.hpp"
19#include "ClpLinearObjective.hpp"
20#include "ClpDualRowSteepest.hpp"
21#include "ClpDualRowDantzig.hpp"
22#include "ClpPrimalColumnSteepest.hpp"
23#include "ClpPrimalColumnDantzig.hpp"
24#include "ClpParameters.hpp"
25#include "ClpNetworkMatrix.hpp"
26#include "ClpPlusMinusOneMatrix.hpp"
27#include "MyMessageHandler.hpp"
28
29#include "ClpPresolve.hpp"
30#ifdef CLP_IDIOT
31#include "Idiot.hpp"
32#endif
33
34
35//#############################################################################
36
37#ifdef NDEBUG
38#undef NDEBUG
39#endif
40
41// Function Prototypes. Function definitions is in this file.
42void testingMessage( const char * const msg );
43
44//----------------------------------------------------------------
45// unitTest [-mpsDir=V1] [-netlibDir=V2] [-netlib]
46//
47// where:
48//   -mpsDir: directory containing mps test files
49//       Default value V1="../Mps/Sample"   
50//   -netlibDir: directory containing netlib files
51//       Default value V2="../Mps/Netlib"
52//   -netlib
53//       If specified, then netlib test set run
54//
55// All parameters are optional.
56//----------------------------------------------------------------
57int mainTest (int argc, const char *argv[],bool doDual,
58              ClpSimplex empty, bool doPresolve,int doIdiot)
59{
60  int i;
61
62  // define valid parameter keywords
63  std::set<std::string> definedKeyWords;
64  definedKeyWords.insert("-mpsDir");
65  definedKeyWords.insert("-netlibDir");
66  definedKeyWords.insert("-netlib");
67
68  // Create a map of parameter keys and associated data
69  std::map<std::string,std::string> parms;
70  for ( i=1; i<argc; i++ ) {
71    std::string parm(argv[i]);
72    std::string key,value;
73    unsigned int  eqPos = parm.find('=');
74
75    // Does parm contain and '='
76    if ( eqPos==std::string::npos ) {
77      //Parm does not contain '='
78      key = parm;
79    }
80    else {
81      key=parm.substr(0,eqPos);
82      value=parm.substr(eqPos+1);
83    }
84
85    // Is specifed key valid?
86    if ( definedKeyWords.find(key) == definedKeyWords.end() ) {
87      // invalid key word.
88      // Write help text
89      std::cerr <<"Undefined parameter \"" <<key <<"\".\n";
90      std::cerr <<"Correct usage: \n";
91      std::cerr <<"  unitTest [-mpsDir=V1] [-netlibDir=V2] [-netlib]\n";
92      std::cerr <<"  where:\n";
93      std::cerr <<"    -mpsDir: directory containing mps test files\n";
94      std::cerr <<"        Default value V1=\"../Mps/Sample\"\n";
95      std::cerr <<"    -netlibDir: directory containing netlib files\n";
96      std::cerr <<"        Default value V2=\"../Mps/Netlib\"\n";
97      std::cerr <<"    -netlib\n";
98      std::cerr <<"        If specified, then netlib testset run.\n";
99      return 1;
100    }
101    parms[key]=value;
102  }
103 
104  const char dirsep =  CoinFindDirSeparator();
105  // Set directory containing mps data files.
106  std::string mpsDir;
107  if (parms.find("-mpsDir") != parms.end())
108    mpsDir=parms["-mpsDir"] + dirsep;
109  else 
110    mpsDir = dirsep == '/' ? "../Mps/Sample/" : "..\\Mps\\Sample\\";
111 
112  // Set directory containing netlib data files.
113  std::string netlibDir;
114  if (parms.find("-netlibDir") != parms.end())
115    netlibDir=parms["-netlibDir"] + dirsep;
116  else 
117    netlibDir = dirsep == '/' ? "../Mps/Netlib/" : "..\\Mps\\Netlib\\";
118
119  testingMessage( "Testing ClpSimplex\n" );
120  ClpSimplexUnitTest(mpsDir,netlibDir);
121  if (parms.find("-netlib") != parms.end())
122  {
123    unsigned int m;
124   
125    // Define test problems:
126    //   mps names,
127    //   maximization or minimization,
128    //   Number of rows and columns in problem, and
129    //   objective function value
130    std::vector<std::string> mpsName;
131    std::vector<bool> min;
132    std::vector<int> nRows;
133    std::vector<int> nCols;
134    std::vector<double> objValue;
135    std::vector<double> objValueTol;
136    mpsName.push_back("25fv47");
137    min.push_back(true);
138    nRows.push_back(822);
139    nCols.push_back(1571);
140    objValueTol.push_back(1.E-10);
141    objValue.push_back(5.5018458883E+03);
142    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);
143    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);
144    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);
145    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);
146   
147    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);
148    mpsName.push_back("pilot4");min.push_back(true);nRows.push_back(411);nCols.push_back(1000);objValueTol.push_back(1.e-8);objValue.push_back(-2.5811392641e+03);
149    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);
150    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);
151    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);
152    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);
153    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);
154    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);
155    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);
156    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);
157    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);
158    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);
159    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);
160    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);
161    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);
162    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);
163    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);
164    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);
165    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);
166    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);
167    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);
168    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);
169    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);
170    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);
171    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); // The correct answer includes -7.113 term. This is a constant in the objective function. See line 1683 of the mps file.
172    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 );
173    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);
174    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);
175    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);
176    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);
177    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);
178    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);
179    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);
180    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);
181    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);
182    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);
183    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);
184    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);
185    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);
186    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);
187    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);
188    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);
189    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);
190    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);
191    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);
192    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);
193    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);
194    //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);
195    //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);
196    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);
197    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);
198    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);
199    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);
200    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);
201    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);
202    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);
203    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);
204    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);
205    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);
206    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);
207    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);
208    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);
209    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);
210    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);
211    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);
212    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);
213    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);
214    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);
215    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);
216    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);
217    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);
218    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);
219    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);
220    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);
221    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);
222    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);
223    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);
224    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);
225    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);
226    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);
227    //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);
228    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); 
229    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);
230    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);
231    //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);
232    //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);
233    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);
234    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);
235    mpsName.push_back("wood1p");min.push_back(true);nRows.push_back(245);nCols.push_back(2594);objValueTol.push_back(1.e-10);objValue.push_back(1.4429024116e+00);
236    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);
237
238    double timeTaken =0.0;
239  // Loop once for each Mps File
240    for (m=0; m<mpsName.size(); m++ ) {
241      std::cerr <<"  processing mps file: " <<mpsName[m] 
242                <<" (" <<m+1 <<" out of " <<mpsName.size() <<")" <<std::endl;
243   
244      // Read data mps file,
245      std::string fn = netlibDir+mpsName[m];
246      CoinMpsIO mps;
247      mps.readMps(fn.c_str(),"mps");
248      double time1 = CoinCpuTime();
249      ClpSimplex solution=empty;
250      solution.loadProblem(*mps.getMatrixByCol(),mps.getColLower(),
251                           mps.getColUpper(),
252                           mps.getObjCoefficients(),
253                           mps.getRowLower(),mps.getRowUpper());
254
255      solution.setDblParam(ClpObjOffset,mps.objectiveOffset());
256#if 0
257      solution.setOptimizationDirection(-1);
258      {
259        int j;
260        double * obj = solution.objective();
261        int n=solution.numberColumns();
262        for (j=0;j<n;j++)
263          obj[j] *= -1.0;
264      }
265#endif
266      if (doPresolve) {
267#ifdef USE_PRESOLVE
268        ClpPresolve pinfo;
269        double a=CoinCpuTime();
270        ClpSimplex * model2 = pinfo.presolvedModel(solution,1.0e-8,
271                                                   false,5,true);
272        printf("Presolve took %g seconds\n",CoinCpuTime()-a);
273        // change from 200 (unless user has changed)
274        if (model2->factorization()->maximumPivots()==200&&
275            model2->numberRows()>-5000)
276          model2->factorization()->maximumPivots(100+model2->numberRows()/100);
277        if (doDual) {
278          // faster if bounds tightened
279          int numberInfeasibilities = model2->tightenPrimalBounds();
280          if (numberInfeasibilities)
281            std::cout<<"** Analysis indicates model infeasible"
282                     <<std::endl;
283          if (doIdiot<0)
284            model2->crash(1000,1);
285          model2->dual();
286        } else {
287#ifdef CLP_IDIOT
288          if (doIdiot>0) {
289            Idiot info(*model2);
290            info.crash(doIdiot);
291          }
292#endif
293          model2->primal(1);
294        }
295        pinfo.postsolve(true);
296
297        delete model2;
298        solution.checkSolution();
299        if (!solution.isProvenOptimal()) {
300          printf("Resolving from postsolved model\n");
301         
302          int savePerturbation = solution.perturbation();
303          solution.setPerturbation(100);
304          solution.primal(1);
305          solution.setPerturbation(savePerturbation);
306          if (solution.numberIterations())
307            printf("****** iterated %d\n",solution.numberIterations());
308        }
309        /*solution.checkSolution();
310        printf("%g dual %g(%d) Primal %g(%d)\n",
311               solution.objectiveValue(),
312               solution.sumDualInfeasibilities(),
313               solution.numberDualInfeasibilities(),
314               solution.sumPrimalInfeasibilities(),
315               solution.numberPrimalInfeasibilities());*/
316        if (0) {
317          ClpPresolve pinfoA;
318          model2 = pinfoA.presolvedModel(solution,1.0e-8);
319
320          printf("Resolving from presolved optimal solution\n");
321          model2->primal(1);
322
323          delete model2;
324        }
325#else
326        if (doDual) {
327          if (doIdiot<0)
328            solution.crash(1000,1);
329          solution.dual();
330        } else {
331#ifdef CLP_IDIOT
332          if (doIdiot>0) {
333            Idiot info(solution);
334            info.crash(doIdiot);
335          }
336#endif
337          solution.primal(1);
338        }
339#endif
340      } else {
341        if (doDual) {
342          if (doIdiot<0)
343            solution.crash(1000,1);
344          solution.dual();
345        } else {
346#ifdef CLP_IDIOT
347          if (doIdiot>0) {
348            Idiot info(solution);
349            info.crash(doIdiot);
350          }
351#endif
352          solution.primal(1);
353        }
354      }
355      double time2 = CoinCpuTime()-time1;
356      timeTaken += time2;
357      printf("Took %g seconds\n",time2);
358      // Test objective solution value
359      {
360        double soln = solution.objectiveValue();
361        CoinRelFltEq eq(objValueTol[m]);
362        std::cerr <<soln <<",  " <<objValue[m] <<" diff "<<
363          soln-objValue[m]<<std::endl;
364        if(!eq(soln,objValue[m]))
365          printf("** difference fails\n");
366      }
367    }
368    printf("Total time %g seconds\n",timeTaken);
369  }
370  else {
371    testingMessage( "***Skipped Testing on netlib    ***\n" );
372    testingMessage( "***use -netlib to test class***\n" );
373  }
374 
375  testingMessage( "All tests completed successfully\n" );
376  return 0;
377}
378
379 
380// Display message on stdout and stderr
381void testingMessage( const char * const msg )
382{
383  std::cerr <<msg;
384  //cout <<endl <<"*****************************************"
385  //     <<endl <<msg <<endl;
386}
387
388//--------------------------------------------------------------------------
389// test factorization methods and simplex method
390void
391ClpSimplexUnitTest(const std::string & mpsDir,
392                   const std::string & netlibDir)
393{
394 
395  CoinRelFltEq eq(0.000001);
396
397  {
398    ClpSimplex solution;
399 
400    // matrix data
401    //deliberate hiccup of 2 between 0 and 1
402    CoinBigIndex start[5]={0,4,7,8,9};
403    int length[5]={2,3,1,1,1};
404    int rows[11]={0,2,-1,-1,0,1,2,0,1,2};
405    double elements[11]={7.0,2.0,1.0e10,1.0e10,-2.0,1.0,-2.0,1,1,1};
406    CoinPackedMatrix matrix(true,3,5,8,elements,rows,start,length);
407   
408    // rim data
409    double objective[7]={-4.0,1.0,0.0,0.0,0.0,0.0,0.0};
410    double rowLower[5]={14.0,3.0,3.0,1.0e10,1.0e10};
411    double rowUpper[5]={14.0,3.0,3.0,-1.0e10,-1.0e10};
412    double colLower[7]={0.0,0.0,0.0,0.0,0.0,0.0,0.0};
413    double colUpper[7]={100.0,100.0,100.0,100.0,100.0,100.0,100.0};
414   
415    // basis 1
416    int rowBasis1[3]={-1,-1,-1};
417    int colBasis1[5]={1,1,-1,-1,1};
418    solution.loadProblem(matrix,colLower,colUpper,objective,
419                         rowLower,rowUpper);
420    int i;
421    solution.createStatus();
422    for (i=0;i<3;i++) {
423      if (rowBasis1[i]<0) {
424        solution.setRowStatus(i,ClpSimplex::atLowerBound);
425      } else {
426        solution.setRowStatus(i,ClpSimplex::basic);
427      }
428    }
429    for (i=0;i<5;i++) {
430      if (colBasis1[i]<0) {
431        solution.setColumnStatus(i,ClpSimplex::atLowerBound);
432      } else {
433        solution.setColumnStatus(i,ClpSimplex::basic);
434      }
435    }
436    solution.setLogLevel(3+4+8+16+32);
437    solution.primal();
438    for (i=0;i<3;i++) {
439      if (rowBasis1[i]<0) {
440        solution.setRowStatus(i,ClpSimplex::atLowerBound);
441      } else {
442        solution.setRowStatus(i,ClpSimplex::basic);
443      }
444    }
445    for (i=0;i<5;i++) {
446      if (colBasis1[i]<0) {
447        solution.setColumnStatus(i,ClpSimplex::atLowerBound);
448      } else {
449        solution.setColumnStatus(i,ClpSimplex::basic);
450      }
451    }
452    // intricate stuff does not work with scaling
453    solution.scaling(0);
454    assert(!solution.factorize ( ));
455    const double * colsol = solution.primalColumnSolution();
456    const double * rowsol = solution.primalRowSolution();
457    solution.getSolution(rowsol,colsol);
458    double colsol1[5]={20.0/7.0,3.0,0.0,0.0,23.0/7.0};
459    for (i=0;i<5;i++) {
460      assert(eq(colsol[i],colsol1[i]));
461    }
462    // now feed in again without actually doing factorization
463    ClpFactorization factorization2 = *solution.factorization();
464    ClpSimplex solution2 = solution;
465    solution2.setFactorization(factorization2);
466    solution2.createStatus();
467    for (i=0;i<3;i++) {
468      if (rowBasis1[i]<0) {
469        solution2.setRowStatus(i,ClpSimplex::atLowerBound);
470      } else {
471        solution2.setRowStatus(i,ClpSimplex::basic);
472      }
473    }
474    for (i=0;i<5;i++) {
475      if (colBasis1[i]<0) {
476        solution2.setColumnStatus(i,ClpSimplex::atLowerBound);
477      } else {
478        solution2.setColumnStatus(i,ClpSimplex::basic);
479      }
480    }
481    // intricate stuff does not work with scaling
482    solution2.scaling(0);
483    solution2.getSolution(rowsol,colsol);
484    colsol = solution2.primalColumnSolution();
485    rowsol = solution2.primalRowSolution();
486    for (i=0;i<5;i++) {
487      assert(eq(colsol[i],colsol1[i]));
488    }
489    solution2.setDualBound(0.1);
490    solution2.dual();
491    objective[2]=-1.0;
492    objective[3]=-0.5;
493    objective[4]=10.0;
494    solution.dual();
495    for (i=0;i<3;i++) {
496      rowLower[i]=-1.0e50;
497      colUpper[i+2]=0.0;
498    }
499    solution.setLogLevel(3);
500    solution.dual();
501    double rowObjective[]={1.0,0.5,-10.0};
502    solution.loadProblem(matrix,colLower,colUpper,objective,
503                         rowLower,rowUpper,rowObjective);
504    solution.dual();
505  }
506  {   
507    CoinMpsIO m;
508    std::string fn = mpsDir+"exmip1";
509    m.readMps(fn.c_str(),"mps");
510    ClpSimplex solution;
511    solution.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
512                         m.getObjCoefficients(),
513                         m.getRowLower(),m.getRowUpper());
514    solution.dual();
515  }
516  // Test Message handler
517  {   
518    CoinMpsIO m;
519    std::string fn = mpsDir+"exmip1";
520    m.readMps(fn.c_str(),"mps");
521    ClpSimplex solution;
522    solution.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
523                         m.getObjCoefficients(),
524                         m.getRowLower(),m.getRowUpper());
525    // Message handler
526    MyMessageHandler messageHandler(&solution);
527    std::cout<<"Testing derived message handler"<<std::endl;
528    solution.passInMessageHandler(&messageHandler);
529    solution.primal();
530
531    // Write saved solutions
532    int nc = solution.getNumCols();
533    int s; 
534    std::deque<StdVectorDouble> fep = messageHandler.getFeasibleExtremePoints();
535    int numSavedSolutions = fep.size();
536    for ( s=0; s<numSavedSolutions; ++s ) {
537      const StdVectorDouble & solnVec = fep[s];
538      for ( int c=0; c<nc; ++c ) {
539        if (fabs(solnVec[c])>1.0e-8)
540          std::cout <<"Saved Solution: " <<s <<" ColNum: " <<c <<" Value: " <<solnVec[c] <<std::endl;
541      }
542    }
543    // solve again without scaling
544    messageHandler.clearFeasibleExtremePoints();
545    solution.scaling(0);
546    solution.allSlackBasis();
547    solution.primal();
548    fep = messageHandler.getFeasibleExtremePoints();
549    numSavedSolutions = fep.size();
550    for ( s=0; s<numSavedSolutions; ++s ) {
551      const StdVectorDouble & solnVec = fep[s];
552      for ( int c=0; c<nc; ++c ) {
553        if (fabs(solnVec[c])>1.0e-8)
554          std::cout <<"Saved Solution: " <<s <<" ColNum: " <<c <<" Value: " <<solnVec[c] <<std::endl;
555      }
556    }
557  }
558  // test steepest edge
559  {   
560    CoinMpsIO m;
561    std::string fn = netlibDir+"finnis";
562    m.readMps(fn.c_str(),"mps");
563    ClpModel model;
564    model.loadProblem(*m.getMatrixByCol(),m.getColLower(),
565                    m.getColUpper(),
566                    m.getObjCoefficients(),
567                    m.getRowLower(),m.getRowUpper());
568    ClpSimplex solution(model);
569
570    solution.scaling(1); 
571    solution.setDualBound(1.0e8);
572    //solution.factorization()->maximumPivots(1);
573    //solution.setLogLevel(3);
574    solution.setDualTolerance(1.0e-7);
575    // set objective sense,
576    ClpDualRowSteepest steep;
577    solution.setDualRowPivotAlgorithm(steep);
578    solution.setDblParam(ClpObjOffset,m.objectiveOffset());
579    solution.dual();
580  }
581  // test normal solution
582  {   
583    CoinMpsIO m;
584    std::string fn = netlibDir+"afiro";
585    m.readMps(fn.c_str(),"mps");
586    ClpSimplex solution;
587    ClpModel model;
588    // do twice - without and with scaling
589    int iPass;
590    for (iPass=0;iPass<2;iPass++) {
591      // explicit row objective for testing
592      int nr = m.getNumRows();
593      double * rowObj = new double[nr];
594      CoinFillN(rowObj,nr,0.0);
595      model.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
596                      m.getObjCoefficients(),
597                      m.getRowLower(),m.getRowUpper(),rowObj);
598      delete [] rowObj;
599      solution = ClpSimplex(model);
600      if (iPass) {
601        solution.scaling();
602      }
603      solution.dual();
604      solution.dual();
605      // test optimal
606      assert (solution.status()==0);
607      int numberColumns = solution.numberColumns();
608      int numberRows = solution.numberRows();
609      CoinPackedVector colsol(numberColumns,solution.primalColumnSolution());
610      double * objective = solution.objective();
611      double objValue = colsol.dotProduct(objective);
612      CoinRelFltEq eq(1.0e-8);
613      assert(eq(objValue,-4.6475314286e+02));
614      double * lower = solution.columnLower();
615      double * upper = solution.columnUpper();
616      double * sol = solution.primalColumnSolution();
617      double * result = new double[numberColumns];
618      CoinFillN ( result, numberColumns,0.0);
619      solution.matrix()->transposeTimes(solution.dualRowSolution(), result);
620      int iRow , iColumn;
621      // see if feasible and dual feasible
622      for (iColumn=0;iColumn<numberColumns;iColumn++) {
623        double value = sol[iColumn];
624        assert(value<upper[iColumn]+1.0e-8);
625        assert(value>lower[iColumn]-1.0e-8);
626        value = objective[iColumn]-result[iColumn];
627        assert (value>-1.0e-5);
628        if (sol[iColumn]>1.0e-5)
629          assert (value<1.0e-5);
630      }
631      delete [] result;
632      result = new double[numberRows];
633      CoinFillN ( result, numberRows,0.0);
634      solution.matrix()->times(colsol, result);
635      lower = solution.rowLower();
636      upper = solution.rowUpper();
637      sol = solution.primalRowSolution();
638      for (iRow=0;iRow<numberRows;iRow++) {
639        double value = result[iRow];
640        assert(eq(value,sol[iRow]));
641        assert(value<upper[iRow]+1.0e-8);
642        assert(value>lower[iRow]-1.0e-8);
643      }
644      delete [] result;
645      // test row objective
646      double * rowObjective = solution.rowObjective();
647      CoinDisjointCopyN(solution.dualRowSolution(),numberRows,rowObjective);
648      CoinDisjointCopyN(solution.dualColumnSolution(),numberColumns,objective);
649      // this sets up all slack basis
650      solution.createStatus();
651      solution.dual();
652      CoinFillN(rowObjective,numberRows,0.0);
653      CoinDisjointCopyN(m.getObjCoefficients(),numberColumns,objective);
654      solution.dual();
655    }
656  }
657  // test unbounded
658  {   
659    CoinMpsIO m;
660    std::string fn = netlibDir+"brandy";
661    m.readMps(fn.c_str(),"mps");
662    ClpSimplex solution;
663    // do twice - without and with scaling
664    int iPass;
665    for (iPass=0;iPass<2;iPass++) {
666      solution.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
667                      m.getObjCoefficients(),
668                      m.getRowLower(),m.getRowUpper());
669      if (iPass)
670        solution.scaling();
671      solution.setOptimizationDirection(-1);
672      // test unbounded and ray
673#ifdef DUAL
674      solution.setDualBound(100.0);
675      solution.dual();
676#else
677      solution.primal();
678#endif
679      assert (solution.status()==2);
680      int numberColumns = solution.numberColumns();
681      int numberRows = solution.numberRows();
682      double * lower = solution.columnLower();
683      double * upper = solution.columnUpper();
684      double * sol = solution.primalColumnSolution();
685      double * ray = solution.unboundedRay();
686      double * objective = solution.objective();
687      double objChange=0.0;
688      int iRow , iColumn;
689      // make sure feasible and columns form ray
690      for (iColumn=0;iColumn<numberColumns;iColumn++) {
691        double value = sol[iColumn];
692        assert(value<upper[iColumn]+1.0e-8);
693        assert(value>lower[iColumn]-1.0e-8);
694        value = ray[iColumn];
695        if (value>0.0)
696          assert(upper[iColumn]>1.0e30);
697        else if (value<0.0)
698          assert(lower[iColumn]<-1.0e30);
699        objChange += value*objective[iColumn];
700      }
701      // make sure increasing objective
702      assert(objChange>0.0);
703      double * result = new double[numberRows];
704      CoinFillN ( result, numberRows,0.0);
705      solution.matrix()->times(sol, result);
706      lower = solution.rowLower();
707      upper = solution.rowUpper();
708      sol = solution.primalRowSolution();
709      for (iRow=0;iRow<numberRows;iRow++) {
710        double value = result[iRow];
711        assert(eq(value,sol[iRow]));
712        assert(value<upper[iRow]+1.0e-8);
713        assert(value>lower[iRow]-1.0e-8);
714      }
715      CoinFillN ( result, numberRows,0.0);
716      solution.matrix()->times(ray, result);
717      // there may be small differences (especially if scaled)
718      for (iRow=0;iRow<numberRows;iRow++) {
719        double value = result[iRow];
720        if (value>1.0e-8)
721          assert(upper[iRow]>1.0e30);
722        else if (value<-1.0e-8)
723          assert(lower[iRow]<-1.0e30);
724      }
725      delete [] result;
726      delete [] ray;
727    }
728  }
729  // test infeasible
730  {   
731    CoinMpsIO m;
732    std::string fn = netlibDir+"brandy";
733    m.readMps(fn.c_str(),"mps");
734    ClpSimplex solution;
735    // do twice - without and with scaling
736    int iPass;
737    for (iPass=0;iPass<2;iPass++) {
738      solution.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
739                      m.getObjCoefficients(),
740                      m.getRowLower(),m.getRowUpper());
741      if (iPass)
742        solution.scaling();
743      // test infeasible and ray
744      solution.columnUpper()[0]=0.0;
745#ifdef DUAL
746      solution.setDualBound(100.0);
747      solution.dual();
748#else
749      solution.primal();
750#endif
751      assert (solution.status()==1);
752      int numberColumns = solution.numberColumns();
753      int numberRows = solution.numberRows();
754      double * lower = solution.rowLower();
755      double * upper = solution.rowUpper();
756      double * ray = solution.infeasibilityRay();
757      assert(ray);
758      // construct proof of infeasibility
759      int iRow , iColumn;
760      double lo=0.0,up=0.0;
761      int nl=0,nu=0;
762      for (iRow=0;iRow<numberRows;iRow++) {
763        if (lower[iRow]>-1.0e20) {
764          lo += ray[iRow]*lower[iRow];
765        } else {
766          if (ray[iRow]>1.0e-8) 
767            nl++;
768        }
769        if (upper[iRow]<1.0e20) {
770          up += ray[iRow]*upper[iRow];
771        } else {
772          if (ray[iRow]>1.0e-8) 
773            nu++;
774        }
775      }
776      if (nl)
777        lo=-1.0e100;
778      if (nu)
779        up=1.0e100;
780      double * result = new double[numberColumns];
781      double lo2=0.0,up2=0.0;
782      CoinFillN ( result, numberColumns,0.0);
783      solution.matrix()->transposeTimes(ray, result);
784      lower = solution.columnLower();
785      upper = solution.columnUpper();
786      nl=nu=0;
787      for (iColumn=0;iColumn<numberColumns;iColumn++) {
788        if (result[iColumn]>1.0e-8) {
789          if (lower[iColumn]>-1.0e20)
790            lo2 += result[iColumn]*lower[iColumn];
791          else
792            nl++;
793          if (upper[iColumn]<1.0e20)
794            up2 += result[iColumn]*upper[iColumn];
795          else
796            nu++;
797        } else if (result[iColumn]<-1.0e-8) {
798          if (lower[iColumn]>-1.0e20)
799            up2 += result[iColumn]*lower[iColumn];
800          else
801            nu++;
802          if (upper[iColumn]<1.0e20)
803            lo2 += result[iColumn]*upper[iColumn];
804          else
805            nl++;
806        }
807      }
808      if (nl)
809        lo2=-1.0e100;
810      if (nu)
811        up2=1.0e100;
812      // make sure inconsistency
813      assert(lo2>up||up2<lo);
814      delete [] result;
815      delete [] ray;
816    }
817  }
818  // test delete and add
819  {   
820    CoinMpsIO m;
821    std::string fn = netlibDir+"brandy";
822    m.readMps(fn.c_str(),"mps");
823    ClpSimplex solution;
824    solution.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
825                         m.getObjCoefficients(),
826                         m.getRowLower(),m.getRowUpper());
827    solution.dual();
828    CoinRelFltEq eq(1.0e-8);
829    assert(eq(solution.objectiveValue(),1.5185098965e+03));
830
831    int numberColumns = solution.numberColumns();
832    int numberRows = solution.numberRows();
833    double * saveObj = new double [numberColumns];
834    double * saveLower = new double[numberRows+numberColumns];
835    double * saveUpper = new double[numberRows+numberColumns];
836    int * which = new int [numberRows+numberColumns];
837
838    int numberElements = m.getMatrixByCol()->getNumElements();
839    int * starts = new int[numberRows+numberColumns];
840    int * index = new int[numberElements];
841    double * element = new double[numberElements];
842
843    const int * startM;
844    const int * lengthM;
845    const int * indexM;
846    const double * elementM;
847
848    int n,nel;
849
850    // delete non basic columns
851    n=0;
852    nel=0;
853    int iRow , iColumn;
854    double * lower = solution.columnLower();
855    double * upper = solution.columnUpper();
856    double * objective = solution.objective();
857    startM = m.getMatrixByCol()->getVectorStarts();
858    lengthM = m.getMatrixByCol()->getVectorLengths();
859    indexM = m.getMatrixByCol()->getIndices();
860    elementM = m.getMatrixByCol()->getElements();
861    starts[0]=0;
862    for (iColumn=0;iColumn<numberColumns;iColumn++) {
863      if (solution.getColumnStatus(iColumn)!=ClpSimplex::basic) {
864        saveObj[n]=objective[iColumn];
865        saveLower[n]=lower[iColumn];
866        saveUpper[n]=upper[iColumn];
867        int j;
868        for (j=startM[iColumn];j<startM[iColumn]+lengthM[iColumn];j++) {
869          index[nel]=indexM[j];
870          element[nel++]=elementM[j];
871        }
872        which[n++]=iColumn;
873        starts[n]=nel;
874      }
875    }
876    solution.deleteColumns(n,which);
877    solution.dual();
878    // Put back
879    solution.addColumns(n,saveLower,saveUpper,saveObj,
880                        starts,index,element);
881    solution.dual();
882    assert(eq(solution.objectiveValue(),1.5185098965e+03));
883
884    // reload with original
885    solution.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
886                         m.getObjCoefficients(),
887                         m.getRowLower(),m.getRowUpper());
888    // delete half rows
889    n=0;
890    nel=0;
891    lower = solution.rowLower();
892    upper = solution.rowUpper();
893    startM = m.getMatrixByRow()->getVectorStarts();
894    lengthM = m.getMatrixByRow()->getVectorLengths();
895    indexM = m.getMatrixByRow()->getIndices();
896    elementM = m.getMatrixByRow()->getElements();
897    starts[0]=0;
898    for (iRow=0;iRow<numberRows;iRow++) {
899      if ((iRow&1)==0) {
900        saveLower[n]=lower[iRow];
901        saveUpper[n]=upper[iRow];
902        int j;
903        for (j=startM[iRow];j<startM[iRow]+lengthM[iRow];j++) {
904          index[nel]=indexM[j];
905          element[nel++]=elementM[j];
906        }
907        which[n++]=iRow;
908        starts[n]=nel;
909      }
910    }
911    solution.deleteRows(n,which);
912    solution.dual();
913    // Put back
914    solution.addRows(n,saveLower,saveUpper,
915                        starts,index,element);
916    solution.dual();
917    assert(eq(solution.objectiveValue(),1.5185098965e+03));
918    // Zero out status array to give some interest
919    memset(solution.statusArray()+numberColumns,0,numberRows);
920    solution.primal(1);
921    assert(eq(solution.objectiveValue(),1.5185098965e+03));
922
923    delete [] saveObj;
924    delete [] saveLower;
925    delete [] saveUpper;
926    delete [] which;
927    delete [] starts;
928    delete [] index;
929    delete [] element;
930  }
931  // test network
932#define QUADRATIC
933#ifndef QUADRATIC
934  if (1) {   
935    std::string fn = mpsDir+"input.130";
936    int numberColumns;
937    int numberRows;
938   
939    FILE * fp = fopen(fn.c_str(),"r");
940    if (!fp) {
941      // Try in Samples
942      fn = "Samples/input.130";
943      fp = fopen(fn.c_str(),"r");
944    }
945    if (!fp) {
946      fprintf(stderr,"Unable to open file input.130 in mpsDir or Samples directory\n");
947    } else {
948      int problem;
949      char temp[100];
950      // read and skip
951      fscanf(fp,"%s",temp);
952      assert (!strcmp(temp,"BEGIN"));
953      fscanf(fp,"%*s %*s %d %d %*s %*s %d %*s",&problem, &numberRows, 
954             &numberColumns);
955      // scan down to SUPPLY
956      while (fgets(temp,100,fp)) {
957        if (!strncmp(temp,"SUPPLY",6))
958          break;
959      }
960      if (strncmp(temp,"SUPPLY",6)) {
961        fprintf(stderr,"Unable to find SUPPLY\n");
962        exit(2);
963      }
964      // get space for rhs
965      double * lower = new double[numberRows];
966      double * upper = new double[numberRows];
967      int i;
968      for (i=0;i<numberRows;i++) {
969        lower[i]=0.0;
970        upper[i]=0.0;
971      }
972      // ***** Remember to convert to C notation
973      while (fgets(temp,100,fp)) {
974        int row;
975        int value;
976        if (!strncmp(temp,"ARCS",4))
977          break;
978        sscanf(temp,"%d %d",&row,&value);
979        upper[row-1]=-value;
980        lower[row-1]=-value;
981      }
982      if (strncmp(temp,"ARCS",4)) {
983        fprintf(stderr,"Unable to find ARCS\n");
984        exit(2);
985      }
986      // number of columns may be underestimate
987      int * head = new int[2*numberColumns];
988      int * tail = new int[2*numberColumns];
989      int * ub = new int[2*numberColumns];
990      int * cost = new int[2*numberColumns];
991      // ***** Remember to convert to C notation
992      numberColumns=0;
993      while (fgets(temp,100,fp)) {
994        int iHead;
995        int iTail;
996        int iUb;
997        int iCost;
998        if (!strncmp(temp,"DEMAND",6))
999          break;
1000        sscanf(temp,"%d %d %d %d",&iHead,&iTail,&iCost,&iUb);
1001        iHead--;
1002        iTail--;
1003        head[numberColumns]=iHead;
1004        tail[numberColumns]=iTail;
1005        ub[numberColumns]=iUb;
1006        cost[numberColumns]=iCost;
1007        numberColumns++;
1008      }
1009      if (strncmp(temp,"DEMAND",6)) {
1010        fprintf(stderr,"Unable to find DEMAND\n");
1011        exit(2);
1012      }
1013      // ***** Remember to convert to C notation
1014      while (fgets(temp,100,fp)) {
1015        int row;
1016        int value;
1017        if (!strncmp(temp,"END",3))
1018          break;
1019        sscanf(temp,"%d %d",&row,&value);
1020        upper[row-1]=value;
1021        lower[row-1]=value;
1022      }
1023      if (strncmp(temp,"END",3)) {
1024        fprintf(stderr,"Unable to find END\n");
1025        exit(2);
1026      }
1027      printf("Problem %d has %d rows and %d columns\n",problem,
1028             numberRows,numberColumns);
1029      fclose(fp);
1030      ClpSimplex  model;
1031      // now build model
1032     
1033      double * objective =new double[numberColumns];
1034      double * lowerColumn = new double[numberColumns];
1035      double * upperColumn = new double[numberColumns];
1036     
1037      double * element = new double [2*numberColumns];
1038      int * start = new int[numberColumns+1];
1039      int * row = new int[2*numberColumns];
1040      start[numberColumns]=2*numberColumns;
1041      for (i=0;i<numberColumns;i++) {
1042        start[i]=2*i;
1043        element[2*i]=-1.0;
1044        element[2*i+1]=1.0;
1045        row[2*i]=head[i];
1046        row[2*i+1]=tail[i];
1047        lowerColumn[i]=0.0;
1048        upperColumn[i]=ub[i];
1049        objective[i]=cost[i];
1050      }
1051      // Create Packed Matrix
1052      CoinPackedMatrix matrix;
1053      int * lengths = NULL;
1054      matrix.assignMatrix(true,numberRows,numberColumns,
1055                          2*numberColumns,element,row,start,lengths);
1056      // load model
1057     
1058      model.loadProblem(matrix,
1059                        lowerColumn,upperColumn,objective,
1060                        lower,upper);
1061     
1062      model.factorization()->maximumPivots(200+model.numberRows()/100);
1063      model.createStatus();
1064      double time1 = CoinCpuTime();
1065      model.dual();
1066      std::cout<<"Network problem, ClpPackedMatrix took "<<CoinCpuTime()-time1<<" seconds"<<std::endl;
1067      ClpPlusMinusOneMatrix plusMinus(matrix);
1068      assert (plusMinus.getIndices()); // would be zero if not +- one
1069      model.loadProblem(plusMinus,
1070                        lowerColumn,upperColumn,objective,
1071                        lower,upper);
1072     
1073      model.factorization()->maximumPivots(200+model.numberRows()/100);
1074      model.createStatus();
1075      time1 = CoinCpuTime();
1076      model.dual();
1077      std::cout<<"Network problem, ClpPlusMinusOneMatrix took "<<CoinCpuTime()-time1<<" seconds"<<std::endl;
1078      ClpNetworkMatrix network(numberColumns,head,tail);
1079      model.loadProblem(network,
1080                        lowerColumn,upperColumn,objective,
1081                        lower,upper);
1082     
1083      model.factorization()->maximumPivots(200+model.numberRows()/100);
1084      model.createStatus();
1085      time1 = CoinCpuTime();
1086      model.dual();
1087      std::cout<<"Network problem, ClpNetworkMatrix took "<<CoinCpuTime()-time1<<" seconds"<<std::endl;
1088      delete [] lower;
1089      delete [] upper;
1090      delete [] head;
1091      delete [] tail;
1092      delete [] ub;
1093      delete [] cost;
1094      delete [] objective;
1095      delete [] lowerColumn;
1096      delete [] upperColumn;
1097    }
1098  }
1099#endif
1100#ifdef QUADRATIC
1101  // Test quadratic to solve linear
1102  if (1) {   
1103    CoinMpsIO m;
1104    std::string fn = mpsDir+"exmip1";
1105    m.readMps(fn.c_str(),"mps");
1106    ClpSimplex solution;
1107    solution.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
1108                         m.getObjCoefficients(),
1109                         m.getRowLower(),m.getRowUpper());
1110    //solution.dual();
1111    // get quadratic part
1112    int numberColumns=solution.numberColumns();
1113    int * start=new int [numberColumns+1];
1114    int * column = new int[numberColumns];
1115    double * element = new double[numberColumns];
1116    int i;
1117    start[0]=0;
1118    for (i=0;i<numberColumns;i++) {
1119      start[i+1]=i+1;
1120      column[i]=i;
1121      element[i]=1.0e-1;
1122      element[i]=0.0;
1123    }
1124    // Load up objective
1125    solution.loadQuadraticObjective(numberColumns,start,column,element);
1126    delete [] start;
1127    delete [] column;
1128    delete [] element;
1129    //solution.quadraticSLP(50,1.0e-4);
1130    double objValue = solution.getObjValue();
1131    CoinRelFltEq eq(1.0e-4);
1132    //assert(eq(objValue,820.0));
1133    solution.setLogLevel(63);
1134    solution.quadraticPrimal(1);
1135    objValue = solution.getObjValue();
1136    assert(eq(objValue,3.2368421));
1137    //exit(77);
1138  }
1139  // Test quadratic
1140  if (1) {   
1141    CoinMpsIO m;
1142    std::string fn = mpsDir+"share2qp";
1143    //fn = "share2qpb";
1144    m.readMps(fn.c_str(),"mps");
1145    ClpSimplex model;
1146    model.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
1147                         m.getObjCoefficients(),
1148                         m.getRowLower(),m.getRowUpper());
1149    model.dual();
1150    // get quadratic part
1151    int * start=NULL;
1152    int * column = NULL;
1153    double * element = NULL;
1154    m.readQuadraticMps(NULL,start,column,element,2);
1155    // Load up objective
1156    model.loadQuadraticObjective(model.numberColumns(),start,column,element);
1157    delete [] start;
1158    delete [] column;
1159    delete [] element;
1160#if 0
1161    model.quadraticSLP(50,1.0e-4);
1162#else
1163    // Get feasible
1164    ClpObjective * saveObjective = model.objectiveAsObject()->clone();
1165    int numberColumns=model.numberColumns();
1166    ClpLinearObjective zeroObjective(NULL,numberColumns);
1167    model.setObjective(&zeroObjective);
1168    model.dual();
1169    model.setObjective(saveObjective);
1170    delete saveObjective;
1171#endif
1172    model.setLogLevel(63);
1173    model.quadraticPrimal(1);
1174    double objValue = model.getObjValue();
1175    CoinRelFltEq eq(1.0e-4);
1176    assert(eq(objValue,-400.92));
1177  }
1178  if (0) {   
1179    CoinMpsIO m;
1180    std::string fn = "./beale";
1181    //fn = "./jensen";
1182    m.readMps(fn.c_str(),"mps");
1183    ClpSimplex solution;
1184    solution.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
1185                         m.getObjCoefficients(),
1186                         m.getRowLower(),m.getRowUpper());
1187    solution.setDblParam(ClpObjOffset,m.objectiveOffset());
1188    solution.dual();
1189    // get quadratic part
1190    int * start=NULL;
1191    int * column = NULL;
1192    double * element = NULL;
1193    m.readQuadraticMps(NULL,start,column,element,2);
1194    // Load up objective
1195    solution.loadQuadraticObjective(solution.numberColumns(),start,column,element);
1196    delete [] start;
1197    delete [] column;
1198    delete [] element;
1199    solution.quadraticPrimal(1);
1200    solution.quadraticSLP(50,1.0e-4);
1201    double objValue = solution.getObjValue();
1202    CoinRelFltEq eq(1.0e-4);
1203    assert(eq(objValue,0.5));
1204    solution.quadraticPrimal();
1205    objValue = solution.getObjValue();
1206    assert(eq(objValue,0.5));
1207  }
1208#endif 
1209}
Note: See TracBrowser for help on using the repository browser.