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

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

For Yiming

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 48.1 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    fn = "Test/subGams4";
521    m.readMps(fn.c_str(),"mps");
522    ClpSimplex model;
523    model.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
524                         m.getObjCoefficients(),
525                         m.getRowLower(),m.getRowUpper());
526    // Message handler
527    MyMessageHandler messageHandler(&model);
528    std::cout<<"Testing derived message handler"<<std::endl;
529    model.passInMessageHandler(&messageHandler);
530    model.messagesPointer()->setDetailMessage(1,102);
531    model.setFactorizationFrequency(10);
532    model.primal();
533
534    // Write saved solutions
535    int nc = model.getNumCols();
536    int s; 
537    std::deque<StdVectorDouble> fep = messageHandler.getFeasibleExtremePoints();
538    int numSavedSolutions = fep.size();
539    for ( s=0; s<numSavedSolutions; ++s ) {
540      const StdVectorDouble & solnVec = fep[s];
541      for ( int c=0; c<nc; ++c ) {
542        if (fabs(solnVec[c])>1.0e-8)
543          std::cout <<"Saved Solution: " <<s <<" ColNum: " <<c <<" Value: " <<solnVec[c] <<std::endl;
544      }
545    }
546    // Solve again without scaling
547    // and maximize then minimize
548    messageHandler.clearFeasibleExtremePoints();
549    model.scaling(0);
550    model.setOptimizationDirection(-1);
551    model.primal();
552    model.setOptimizationDirection(1);
553    model.primal();
554    fep = messageHandler.getFeasibleExtremePoints();
555    numSavedSolutions = fep.size();
556    for ( s=0; s<numSavedSolutions; ++s ) {
557      const StdVectorDouble & solnVec = fep[s];
558      for ( int c=0; c<nc; ++c ) {
559        if (fabs(solnVec[c])>1.0e-8)
560          std::cout <<"Saved Solution: " <<s <<" ColNum: " <<c <<" Value: " <<solnVec[c] <<std::endl;
561      }
562    }
563  }
564  // test steepest edge
565  {   
566    CoinMpsIO m;
567    std::string fn = netlibDir+"finnis";
568    m.readMps(fn.c_str(),"mps");
569    ClpModel model;
570    model.loadProblem(*m.getMatrixByCol(),m.getColLower(),
571                    m.getColUpper(),
572                    m.getObjCoefficients(),
573                    m.getRowLower(),m.getRowUpper());
574    ClpSimplex solution(model);
575
576    solution.scaling(1); 
577    solution.setDualBound(1.0e8);
578    //solution.factorization()->maximumPivots(1);
579    //solution.setLogLevel(3);
580    solution.setDualTolerance(1.0e-7);
581    // set objective sense,
582    ClpDualRowSteepest steep;
583    solution.setDualRowPivotAlgorithm(steep);
584    solution.setDblParam(ClpObjOffset,m.objectiveOffset());
585    solution.dual();
586  }
587  // test normal solution
588  {   
589    CoinMpsIO m;
590    std::string fn = netlibDir+"afiro";
591    m.readMps(fn.c_str(),"mps");
592    ClpSimplex solution;
593    ClpModel model;
594    // do twice - without and with scaling
595    int iPass;
596    for (iPass=0;iPass<2;iPass++) {
597      // explicit row objective for testing
598      int nr = m.getNumRows();
599      double * rowObj = new double[nr];
600      CoinFillN(rowObj,nr,0.0);
601      model.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
602                      m.getObjCoefficients(),
603                      m.getRowLower(),m.getRowUpper(),rowObj);
604      delete [] rowObj;
605      solution = ClpSimplex(model);
606      if (iPass) {
607        solution.scaling();
608      }
609      solution.dual();
610      solution.dual();
611      // test optimal
612      assert (solution.status()==0);
613      int numberColumns = solution.numberColumns();
614      int numberRows = solution.numberRows();
615      CoinPackedVector colsol(numberColumns,solution.primalColumnSolution());
616      double * objective = solution.objective();
617      double objValue = colsol.dotProduct(objective);
618      CoinRelFltEq eq(1.0e-8);
619      assert(eq(objValue,-4.6475314286e+02));
620      double * lower = solution.columnLower();
621      double * upper = solution.columnUpper();
622      double * sol = solution.primalColumnSolution();
623      double * result = new double[numberColumns];
624      CoinFillN ( result, numberColumns,0.0);
625      solution.matrix()->transposeTimes(solution.dualRowSolution(), result);
626      int iRow , iColumn;
627      // see if feasible and dual feasible
628      for (iColumn=0;iColumn<numberColumns;iColumn++) {
629        double value = sol[iColumn];
630        assert(value<upper[iColumn]+1.0e-8);
631        assert(value>lower[iColumn]-1.0e-8);
632        value = objective[iColumn]-result[iColumn];
633        assert (value>-1.0e-5);
634        if (sol[iColumn]>1.0e-5)
635          assert (value<1.0e-5);
636      }
637      delete [] result;
638      result = new double[numberRows];
639      CoinFillN ( result, numberRows,0.0);
640      solution.matrix()->times(colsol, result);
641      lower = solution.rowLower();
642      upper = solution.rowUpper();
643      sol = solution.primalRowSolution();
644      for (iRow=0;iRow<numberRows;iRow++) {
645        double value = result[iRow];
646        assert(eq(value,sol[iRow]));
647        assert(value<upper[iRow]+1.0e-8);
648        assert(value>lower[iRow]-1.0e-8);
649      }
650      delete [] result;
651      // test row objective
652      double * rowObjective = solution.rowObjective();
653      CoinDisjointCopyN(solution.dualRowSolution(),numberRows,rowObjective);
654      CoinDisjointCopyN(solution.dualColumnSolution(),numberColumns,objective);
655      // this sets up all slack basis
656      solution.createStatus();
657      solution.dual();
658      CoinFillN(rowObjective,numberRows,0.0);
659      CoinDisjointCopyN(m.getObjCoefficients(),numberColumns,objective);
660      solution.dual();
661    }
662  }
663  // test unbounded
664  {   
665    CoinMpsIO m;
666    std::string fn = netlibDir+"brandy";
667    m.readMps(fn.c_str(),"mps");
668    ClpSimplex solution;
669    // do twice - without and with scaling
670    int iPass;
671    for (iPass=0;iPass<2;iPass++) {
672      solution.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
673                      m.getObjCoefficients(),
674                      m.getRowLower(),m.getRowUpper());
675      if (iPass)
676        solution.scaling();
677      solution.setOptimizationDirection(-1);
678      // test unbounded and ray
679#ifdef DUAL
680      solution.setDualBound(100.0);
681      solution.dual();
682#else
683      solution.primal();
684#endif
685      assert (solution.status()==2);
686      int numberColumns = solution.numberColumns();
687      int numberRows = solution.numberRows();
688      double * lower = solution.columnLower();
689      double * upper = solution.columnUpper();
690      double * sol = solution.primalColumnSolution();
691      double * ray = solution.unboundedRay();
692      double * objective = solution.objective();
693      double objChange=0.0;
694      int iRow , iColumn;
695      // make sure feasible and columns form ray
696      for (iColumn=0;iColumn<numberColumns;iColumn++) {
697        double value = sol[iColumn];
698        assert(value<upper[iColumn]+1.0e-8);
699        assert(value>lower[iColumn]-1.0e-8);
700        value = ray[iColumn];
701        if (value>0.0)
702          assert(upper[iColumn]>1.0e30);
703        else if (value<0.0)
704          assert(lower[iColumn]<-1.0e30);
705        objChange += value*objective[iColumn];
706      }
707      // make sure increasing objective
708      assert(objChange>0.0);
709      double * result = new double[numberRows];
710      CoinFillN ( result, numberRows,0.0);
711      solution.matrix()->times(sol, result);
712      lower = solution.rowLower();
713      upper = solution.rowUpper();
714      sol = solution.primalRowSolution();
715      for (iRow=0;iRow<numberRows;iRow++) {
716        double value = result[iRow];
717        assert(eq(value,sol[iRow]));
718        assert(value<upper[iRow]+1.0e-8);
719        assert(value>lower[iRow]-1.0e-8);
720      }
721      CoinFillN ( result, numberRows,0.0);
722      solution.matrix()->times(ray, result);
723      // there may be small differences (especially if scaled)
724      for (iRow=0;iRow<numberRows;iRow++) {
725        double value = result[iRow];
726        if (value>1.0e-8)
727          assert(upper[iRow]>1.0e30);
728        else if (value<-1.0e-8)
729          assert(lower[iRow]<-1.0e30);
730      }
731      delete [] result;
732      delete [] ray;
733    }
734  }
735  // test infeasible
736  {   
737    CoinMpsIO m;
738    std::string fn = netlibDir+"brandy";
739    m.readMps(fn.c_str(),"mps");
740    ClpSimplex solution;
741    // do twice - without and with scaling
742    int iPass;
743    for (iPass=0;iPass<2;iPass++) {
744      solution.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
745                      m.getObjCoefficients(),
746                      m.getRowLower(),m.getRowUpper());
747      if (iPass)
748        solution.scaling();
749      // test infeasible and ray
750      solution.columnUpper()[0]=0.0;
751#ifdef DUAL
752      solution.setDualBound(100.0);
753      solution.dual();
754#else
755      solution.primal();
756#endif
757      assert (solution.status()==1);
758      int numberColumns = solution.numberColumns();
759      int numberRows = solution.numberRows();
760      double * lower = solution.rowLower();
761      double * upper = solution.rowUpper();
762      double * ray = solution.infeasibilityRay();
763      assert(ray);
764      // construct proof of infeasibility
765      int iRow , iColumn;
766      double lo=0.0,up=0.0;
767      int nl=0,nu=0;
768      for (iRow=0;iRow<numberRows;iRow++) {
769        if (lower[iRow]>-1.0e20) {
770          lo += ray[iRow]*lower[iRow];
771        } else {
772          if (ray[iRow]>1.0e-8) 
773            nl++;
774        }
775        if (upper[iRow]<1.0e20) {
776          up += ray[iRow]*upper[iRow];
777        } else {
778          if (ray[iRow]>1.0e-8) 
779            nu++;
780        }
781      }
782      if (nl)
783        lo=-1.0e100;
784      if (nu)
785        up=1.0e100;
786      double * result = new double[numberColumns];
787      double lo2=0.0,up2=0.0;
788      CoinFillN ( result, numberColumns,0.0);
789      solution.matrix()->transposeTimes(ray, result);
790      lower = solution.columnLower();
791      upper = solution.columnUpper();
792      nl=nu=0;
793      for (iColumn=0;iColumn<numberColumns;iColumn++) {
794        if (result[iColumn]>1.0e-8) {
795          if (lower[iColumn]>-1.0e20)
796            lo2 += result[iColumn]*lower[iColumn];
797          else
798            nl++;
799          if (upper[iColumn]<1.0e20)
800            up2 += result[iColumn]*upper[iColumn];
801          else
802            nu++;
803        } else if (result[iColumn]<-1.0e-8) {
804          if (lower[iColumn]>-1.0e20)
805            up2 += result[iColumn]*lower[iColumn];
806          else
807            nu++;
808          if (upper[iColumn]<1.0e20)
809            lo2 += result[iColumn]*upper[iColumn];
810          else
811            nl++;
812        }
813      }
814      if (nl)
815        lo2=-1.0e100;
816      if (nu)
817        up2=1.0e100;
818      // make sure inconsistency
819      assert(lo2>up||up2<lo);
820      delete [] result;
821      delete [] ray;
822    }
823  }
824  // test delete and add
825  {   
826    CoinMpsIO m;
827    std::string fn = netlibDir+"brandy";
828    m.readMps(fn.c_str(),"mps");
829    ClpSimplex solution;
830    solution.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
831                         m.getObjCoefficients(),
832                         m.getRowLower(),m.getRowUpper());
833    solution.dual();
834    CoinRelFltEq eq(1.0e-8);
835    assert(eq(solution.objectiveValue(),1.5185098965e+03));
836
837    int numberColumns = solution.numberColumns();
838    int numberRows = solution.numberRows();
839    double * saveObj = new double [numberColumns];
840    double * saveLower = new double[numberRows+numberColumns];
841    double * saveUpper = new double[numberRows+numberColumns];
842    int * which = new int [numberRows+numberColumns];
843
844    int numberElements = m.getMatrixByCol()->getNumElements();
845    int * starts = new int[numberRows+numberColumns];
846    int * index = new int[numberElements];
847    double * element = new double[numberElements];
848
849    const int * startM;
850    const int * lengthM;
851    const int * indexM;
852    const double * elementM;
853
854    int n,nel;
855
856    // delete non basic columns
857    n=0;
858    nel=0;
859    int iRow , iColumn;
860    double * lower = solution.columnLower();
861    double * upper = solution.columnUpper();
862    double * objective = solution.objective();
863    startM = m.getMatrixByCol()->getVectorStarts();
864    lengthM = m.getMatrixByCol()->getVectorLengths();
865    indexM = m.getMatrixByCol()->getIndices();
866    elementM = m.getMatrixByCol()->getElements();
867    starts[0]=0;
868    for (iColumn=0;iColumn<numberColumns;iColumn++) {
869      if (solution.getColumnStatus(iColumn)!=ClpSimplex::basic) {
870        saveObj[n]=objective[iColumn];
871        saveLower[n]=lower[iColumn];
872        saveUpper[n]=upper[iColumn];
873        int j;
874        for (j=startM[iColumn];j<startM[iColumn]+lengthM[iColumn];j++) {
875          index[nel]=indexM[j];
876          element[nel++]=elementM[j];
877        }
878        which[n++]=iColumn;
879        starts[n]=nel;
880      }
881    }
882    solution.deleteColumns(n,which);
883    solution.dual();
884    // Put back
885    solution.addColumns(n,saveLower,saveUpper,saveObj,
886                        starts,index,element);
887    solution.dual();
888    assert(eq(solution.objectiveValue(),1.5185098965e+03));
889
890    // reload with original
891    solution.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
892                         m.getObjCoefficients(),
893                         m.getRowLower(),m.getRowUpper());
894    // delete half rows
895    n=0;
896    nel=0;
897    lower = solution.rowLower();
898    upper = solution.rowUpper();
899    startM = m.getMatrixByRow()->getVectorStarts();
900    lengthM = m.getMatrixByRow()->getVectorLengths();
901    indexM = m.getMatrixByRow()->getIndices();
902    elementM = m.getMatrixByRow()->getElements();
903    starts[0]=0;
904    for (iRow=0;iRow<numberRows;iRow++) {
905      if ((iRow&1)==0) {
906        saveLower[n]=lower[iRow];
907        saveUpper[n]=upper[iRow];
908        int j;
909        for (j=startM[iRow];j<startM[iRow]+lengthM[iRow];j++) {
910          index[nel]=indexM[j];
911          element[nel++]=elementM[j];
912        }
913        which[n++]=iRow;
914        starts[n]=nel;
915      }
916    }
917    solution.deleteRows(n,which);
918    solution.dual();
919    // Put back
920    solution.addRows(n,saveLower,saveUpper,
921                        starts,index,element);
922    solution.dual();
923    assert(eq(solution.objectiveValue(),1.5185098965e+03));
924    // Zero out status array to give some interest
925    memset(solution.statusArray()+numberColumns,0,numberRows);
926    solution.primal(1);
927    assert(eq(solution.objectiveValue(),1.5185098965e+03));
928
929    delete [] saveObj;
930    delete [] saveLower;
931    delete [] saveUpper;
932    delete [] which;
933    delete [] starts;
934    delete [] index;
935    delete [] element;
936  }
937  // test network
938#define QUADRATIC
939#ifndef QUADRATIC
940  if (1) {   
941    std::string fn = mpsDir+"input.130";
942    int numberColumns;
943    int numberRows;
944   
945    FILE * fp = fopen(fn.c_str(),"r");
946    if (!fp) {
947      // Try in Samples
948      fn = "Samples/input.130";
949      fp = fopen(fn.c_str(),"r");
950    }
951    if (!fp) {
952      fprintf(stderr,"Unable to open file input.130 in mpsDir or Samples directory\n");
953    } else {
954      int problem;
955      char temp[100];
956      // read and skip
957      fscanf(fp,"%s",temp);
958      assert (!strcmp(temp,"BEGIN"));
959      fscanf(fp,"%*s %*s %d %d %*s %*s %d %*s",&problem, &numberRows, 
960             &numberColumns);
961      // scan down to SUPPLY
962      while (fgets(temp,100,fp)) {
963        if (!strncmp(temp,"SUPPLY",6))
964          break;
965      }
966      if (strncmp(temp,"SUPPLY",6)) {
967        fprintf(stderr,"Unable to find SUPPLY\n");
968        exit(2);
969      }
970      // get space for rhs
971      double * lower = new double[numberRows];
972      double * upper = new double[numberRows];
973      int i;
974      for (i=0;i<numberRows;i++) {
975        lower[i]=0.0;
976        upper[i]=0.0;
977      }
978      // ***** Remember to convert to C notation
979      while (fgets(temp,100,fp)) {
980        int row;
981        int value;
982        if (!strncmp(temp,"ARCS",4))
983          break;
984        sscanf(temp,"%d %d",&row,&value);
985        upper[row-1]=-value;
986        lower[row-1]=-value;
987      }
988      if (strncmp(temp,"ARCS",4)) {
989        fprintf(stderr,"Unable to find ARCS\n");
990        exit(2);
991      }
992      // number of columns may be underestimate
993      int * head = new int[2*numberColumns];
994      int * tail = new int[2*numberColumns];
995      int * ub = new int[2*numberColumns];
996      int * cost = new int[2*numberColumns];
997      // ***** Remember to convert to C notation
998      numberColumns=0;
999      while (fgets(temp,100,fp)) {
1000        int iHead;
1001        int iTail;
1002        int iUb;
1003        int iCost;
1004        if (!strncmp(temp,"DEMAND",6))
1005          break;
1006        sscanf(temp,"%d %d %d %d",&iHead,&iTail,&iCost,&iUb);
1007        iHead--;
1008        iTail--;
1009        head[numberColumns]=iHead;
1010        tail[numberColumns]=iTail;
1011        ub[numberColumns]=iUb;
1012        cost[numberColumns]=iCost;
1013        numberColumns++;
1014      }
1015      if (strncmp(temp,"DEMAND",6)) {
1016        fprintf(stderr,"Unable to find DEMAND\n");
1017        exit(2);
1018      }
1019      // ***** Remember to convert to C notation
1020      while (fgets(temp,100,fp)) {
1021        int row;
1022        int value;
1023        if (!strncmp(temp,"END",3))
1024          break;
1025        sscanf(temp,"%d %d",&row,&value);
1026        upper[row-1]=value;
1027        lower[row-1]=value;
1028      }
1029      if (strncmp(temp,"END",3)) {
1030        fprintf(stderr,"Unable to find END\n");
1031        exit(2);
1032      }
1033      printf("Problem %d has %d rows and %d columns\n",problem,
1034             numberRows,numberColumns);
1035      fclose(fp);
1036      ClpSimplex  model;
1037      // now build model
1038     
1039      double * objective =new double[numberColumns];
1040      double * lowerColumn = new double[numberColumns];
1041      double * upperColumn = new double[numberColumns];
1042     
1043      double * element = new double [2*numberColumns];
1044      int * start = new int[numberColumns+1];
1045      int * row = new int[2*numberColumns];
1046      start[numberColumns]=2*numberColumns;
1047      for (i=0;i<numberColumns;i++) {
1048        start[i]=2*i;
1049        element[2*i]=-1.0;
1050        element[2*i+1]=1.0;
1051        row[2*i]=head[i];
1052        row[2*i+1]=tail[i];
1053        lowerColumn[i]=0.0;
1054        upperColumn[i]=ub[i];
1055        objective[i]=cost[i];
1056      }
1057      // Create Packed Matrix
1058      CoinPackedMatrix matrix;
1059      int * lengths = NULL;
1060      matrix.assignMatrix(true,numberRows,numberColumns,
1061                          2*numberColumns,element,row,start,lengths);
1062      // load model
1063     
1064      model.loadProblem(matrix,
1065                        lowerColumn,upperColumn,objective,
1066                        lower,upper);
1067     
1068      model.factorization()->maximumPivots(200+model.numberRows()/100);
1069      model.createStatus();
1070      double time1 = CoinCpuTime();
1071      model.dual();
1072      std::cout<<"Network problem, ClpPackedMatrix took "<<CoinCpuTime()-time1<<" seconds"<<std::endl;
1073      ClpPlusMinusOneMatrix plusMinus(matrix);
1074      assert (plusMinus.getIndices()); // would be zero if not +- one
1075      model.loadProblem(plusMinus,
1076                        lowerColumn,upperColumn,objective,
1077                        lower,upper);
1078     
1079      model.factorization()->maximumPivots(200+model.numberRows()/100);
1080      model.createStatus();
1081      time1 = CoinCpuTime();
1082      model.dual();
1083      std::cout<<"Network problem, ClpPlusMinusOneMatrix took "<<CoinCpuTime()-time1<<" seconds"<<std::endl;
1084      ClpNetworkMatrix network(numberColumns,head,tail);
1085      model.loadProblem(network,
1086                        lowerColumn,upperColumn,objective,
1087                        lower,upper);
1088     
1089      model.factorization()->maximumPivots(200+model.numberRows()/100);
1090      model.createStatus();
1091      time1 = CoinCpuTime();
1092      model.dual();
1093      std::cout<<"Network problem, ClpNetworkMatrix took "<<CoinCpuTime()-time1<<" seconds"<<std::endl;
1094      delete [] lower;
1095      delete [] upper;
1096      delete [] head;
1097      delete [] tail;
1098      delete [] ub;
1099      delete [] cost;
1100      delete [] objective;
1101      delete [] lowerColumn;
1102      delete [] upperColumn;
1103    }
1104  }
1105#endif
1106#ifdef QUADRATIC
1107  // Test quadratic to solve linear
1108  if (1) {   
1109    CoinMpsIO m;
1110    std::string fn = mpsDir+"exmip1";
1111    m.readMps(fn.c_str(),"mps");
1112    ClpSimplex solution;
1113    solution.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
1114                         m.getObjCoefficients(),
1115                         m.getRowLower(),m.getRowUpper());
1116    //solution.dual();
1117    // get quadratic part
1118    int numberColumns=solution.numberColumns();
1119    int * start=new int [numberColumns+1];
1120    int * column = new int[numberColumns];
1121    double * element = new double[numberColumns];
1122    int i;
1123    start[0]=0;
1124    for (i=0;i<numberColumns;i++) {
1125      start[i+1]=i+1;
1126      column[i]=i;
1127      element[i]=1.0e-1;
1128      element[i]=0.0;
1129    }
1130    // Load up objective
1131    solution.loadQuadraticObjective(numberColumns,start,column,element);
1132    delete [] start;
1133    delete [] column;
1134    delete [] element;
1135    //solution.quadraticSLP(50,1.0e-4);
1136    double objValue = solution.getObjValue();
1137    CoinRelFltEq eq(1.0e-4);
1138    //assert(eq(objValue,820.0));
1139    solution.setLogLevel(63);
1140    solution.quadraticPrimal(1);
1141    objValue = solution.getObjValue();
1142    assert(eq(objValue,3.2368421));
1143    //exit(77);
1144  }
1145  // Test quadratic
1146  if (1) {   
1147    CoinMpsIO m;
1148    std::string fn = mpsDir+"share2qp";
1149    //fn = "share2qpb";
1150    m.readMps(fn.c_str(),"mps");
1151    ClpSimplex model;
1152    model.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
1153                         m.getObjCoefficients(),
1154                         m.getRowLower(),m.getRowUpper());
1155    model.dual();
1156    // get quadratic part
1157    int * start=NULL;
1158    int * column = NULL;
1159    double * element = NULL;
1160    m.readQuadraticMps(NULL,start,column,element,2);
1161    // Load up objective
1162    model.loadQuadraticObjective(model.numberColumns(),start,column,element);
1163    delete [] start;
1164    delete [] column;
1165    delete [] element;
1166#if 0
1167    model.quadraticSLP(50,1.0e-4);
1168#else
1169    // Get feasible
1170    ClpObjective * saveObjective = model.objectiveAsObject()->clone();
1171    int numberColumns=model.numberColumns();
1172    ClpLinearObjective zeroObjective(NULL,numberColumns);
1173    model.setObjective(&zeroObjective);
1174    model.dual();
1175    model.setObjective(saveObjective);
1176    delete saveObjective;
1177#endif
1178    model.setLogLevel(63);
1179    model.quadraticPrimal(1);
1180    double objValue = model.getObjValue();
1181    CoinRelFltEq eq(1.0e-4);
1182    assert(eq(objValue,-400.92));
1183  }
1184  if (0) {   
1185    CoinMpsIO m;
1186    std::string fn = "./beale";
1187    //fn = "./jensen";
1188    m.readMps(fn.c_str(),"mps");
1189    ClpSimplex solution;
1190    solution.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
1191                         m.getObjCoefficients(),
1192                         m.getRowLower(),m.getRowUpper());
1193    solution.setDblParam(ClpObjOffset,m.objectiveOffset());
1194    solution.dual();
1195    // get quadratic part
1196    int * start=NULL;
1197    int * column = NULL;
1198    double * element = NULL;
1199    m.readQuadraticMps(NULL,start,column,element,2);
1200    // Load up objective
1201    solution.loadQuadraticObjective(solution.numberColumns(),start,column,element);
1202    delete [] start;
1203    delete [] column;
1204    delete [] element;
1205    solution.quadraticPrimal(1);
1206    solution.quadraticSLP(50,1.0e-4);
1207    double objValue = solution.getObjValue();
1208    CoinRelFltEq eq(1.0e-4);
1209    assert(eq(objValue,0.5));
1210    solution.quadraticPrimal();
1211    objValue = solution.getObjValue();
1212    assert(eq(objValue,0.5));
1213  }
1214#endif 
1215}
Note: See TracBrowser for help on using the repository browser.