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

Last change on this file since 210 was 210, checked in by forrest, 17 years ago

Trying

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