source: trunk/Test/unitTest.cpp @ 345

Last change on this file since 345 was 345, checked in by jpfasano, 16 years ago

Modified to compile with MS Visual C++ V6

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 51.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 "ClpInterior.hpp"
21#include "ClpLinearObjective.hpp"
22#include "ClpDualRowSteepest.hpp"
23#include "ClpDualRowDantzig.hpp"
24#include "ClpPrimalColumnSteepest.hpp"
25#include "ClpPrimalColumnDantzig.hpp"
26#include "ClpParameters.hpp"
27#include "ClpNetworkMatrix.hpp"
28#include "ClpPlusMinusOneMatrix.hpp"
29#include "MyMessageHandler.hpp"
30#include "MyEventHandler.hpp"
31
32#include "ClpPresolve.hpp"
33#include "Idiot.hpp"
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(5.e-5);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(5.e-5);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        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          if (doIdiot>0) {
288            Idiot info(*model2);
289            info.crash(doIdiot,model2->messageHandler(),model2->messagesPointer());
290          }
291          model2->primal(1);
292        }
293        pinfo.postsolve(true);
294
295        delete model2;
296        solution.checkSolution();
297        if (!solution.isProvenOptimal()) {
298          printf("Resolving from postsolved model\n");
299         
300          int savePerturbation = solution.perturbation();
301          solution.setPerturbation(100);
302          solution.primal(1);
303          solution.setPerturbation(savePerturbation);
304          if (solution.numberIterations())
305            printf("****** iterated %d\n",solution.numberIterations());
306        }
307        /*solution.checkSolution();
308        printf("%g dual %g(%d) Primal %g(%d)\n",
309               solution.objectiveValue(),
310               solution.sumDualInfeasibilities(),
311               solution.numberDualInfeasibilities(),
312               solution.sumPrimalInfeasibilities(),
313               solution.numberPrimalInfeasibilities());*/
314        if (0) {
315          ClpPresolve pinfoA;
316          model2 = pinfoA.presolvedModel(solution,1.0e-8);
317
318          printf("Resolving from presolved optimal solution\n");
319          model2->primal(1);
320
321          delete model2;
322        }
323      } else {
324        if (doDual) {
325          if (doIdiot<0)
326            solution.crash(1000,1);
327          solution.dual();
328        } else {
329          if (doIdiot>0) {
330            Idiot info(solution);
331            info.crash(doIdiot,solution.messageHandler(),solution.messagesPointer());
332          }
333          solution.primal(1);
334        }
335      }
336      double time2 = CoinCpuTime()-time1;
337      timeTaken += time2;
338      printf("Took %g seconds\n",time2);
339      // Test objective solution value
340      {
341        double soln = solution.objectiveValue();
342        CoinRelFltEq eq(objValueTol[m]);
343        std::cerr <<soln <<",  " <<objValue[m] <<" diff "<<
344          soln-objValue[m]<<std::endl;
345        if(!eq(soln,objValue[m]))
346          printf("** difference fails\n");
347      }
348    }
349    printf("Total time %g seconds\n",timeTaken);
350  }
351  else {
352    testingMessage( "***Skipped Testing on netlib    ***\n" );
353    testingMessage( "***use -netlib to test class***\n" );
354  }
355 
356  testingMessage( "All tests completed successfully\n" );
357  return 0;
358}
359
360 
361// Display message on stdout and stderr
362void testingMessage( const char * const msg )
363{
364  std::cerr <<msg;
365  //cout <<endl <<"*****************************************"
366  //     <<endl <<msg <<endl;
367}
368
369//--------------------------------------------------------------------------
370// test factorization methods and simplex method and simple barrier
371void
372ClpSimplexUnitTest(const std::string & mpsDir,
373                   const std::string & netlibDir)
374{
375 
376  CoinRelFltEq eq(0.000001);
377
378  {
379    ClpSimplex solution;
380 
381    // matrix data
382    //deliberate hiccup of 2 between 0 and 1
383    CoinBigIndex start[5]={0,4,7,8,9};
384    int length[5]={2,3,1,1,1};
385    int rows[11]={0,2,-1,-1,0,1,2,0,1,2};
386    double elements[11]={7.0,2.0,1.0e10,1.0e10,-2.0,1.0,-2.0,1,1,1};
387    CoinPackedMatrix matrix(true,3,5,8,elements,rows,start,length);
388   
389    // rim data
390    double objective[7]={-4.0,1.0,0.0,0.0,0.0,0.0,0.0};
391    double rowLower[5]={14.0,3.0,3.0,1.0e10,1.0e10};
392    double rowUpper[5]={14.0,3.0,3.0,-1.0e10,-1.0e10};
393    double colLower[7]={0.0,0.0,0.0,0.0,0.0,0.0,0.0};
394    double colUpper[7]={100.0,100.0,100.0,100.0,100.0,100.0,100.0};
395   
396    // basis 1
397    int rowBasis1[3]={-1,-1,-1};
398    int colBasis1[5]={1,1,-1,-1,1};
399    solution.loadProblem(matrix,colLower,colUpper,objective,
400                         rowLower,rowUpper);
401    int i;
402    solution.createStatus();
403    for (i=0;i<3;i++) {
404      if (rowBasis1[i]<0) {
405        solution.setRowStatus(i,ClpSimplex::atLowerBound);
406      } else {
407        solution.setRowStatus(i,ClpSimplex::basic);
408      }
409    }
410    for (i=0;i<5;i++) {
411      if (colBasis1[i]<0) {
412        solution.setColumnStatus(i,ClpSimplex::atLowerBound);
413      } else {
414        solution.setColumnStatus(i,ClpSimplex::basic);
415      }
416    }
417    solution.setLogLevel(3+4+8+16+32);
418    solution.primal();
419    for (i=0;i<3;i++) {
420      if (rowBasis1[i]<0) {
421        solution.setRowStatus(i,ClpSimplex::atLowerBound);
422      } else {
423        solution.setRowStatus(i,ClpSimplex::basic);
424      }
425    }
426    for (i=0;i<5;i++) {
427      if (colBasis1[i]<0) {
428        solution.setColumnStatus(i,ClpSimplex::atLowerBound);
429      } else {
430        solution.setColumnStatus(i,ClpSimplex::basic);
431      }
432    }
433    // intricate stuff does not work with scaling
434    solution.scaling(0);
435    int returnCode = solution.factorize ( );
436    assert(!returnCode);
437    const double * colsol = solution.primalColumnSolution();
438    const double * rowsol = solution.primalRowSolution();
439    solution.getSolution(rowsol,colsol);
440    double colsol1[5]={20.0/7.0,3.0,0.0,0.0,23.0/7.0};
441    for (i=0;i<5;i++) {
442      assert(eq(colsol[i],colsol1[i]));
443    }
444    // now feed in again without actually doing factorization
445    ClpFactorization factorization2 = *solution.factorization();
446    ClpSimplex solution2 = solution;
447    solution2.setFactorization(factorization2);
448    solution2.createStatus();
449    for (i=0;i<3;i++) {
450      if (rowBasis1[i]<0) {
451        solution2.setRowStatus(i,ClpSimplex::atLowerBound);
452      } else {
453        solution2.setRowStatus(i,ClpSimplex::basic);
454      }
455    }
456    for (i=0;i<5;i++) {
457      if (colBasis1[i]<0) {
458        solution2.setColumnStatus(i,ClpSimplex::atLowerBound);
459      } else {
460        solution2.setColumnStatus(i,ClpSimplex::basic);
461      }
462    }
463    // intricate stuff does not work with scaling
464    solution2.scaling(0);
465    solution2.getSolution(rowsol,colsol);
466    colsol = solution2.primalColumnSolution();
467    rowsol = solution2.primalRowSolution();
468    for (i=0;i<5;i++) {
469      assert(eq(colsol[i],colsol1[i]));
470    }
471    solution2.setDualBound(0.1);
472    solution2.dual();
473    objective[2]=-1.0;
474    objective[3]=-0.5;
475    objective[4]=10.0;
476    solution.dual();
477    for (i=0;i<3;i++) {
478      rowLower[i]=-1.0e50;
479      colUpper[i+2]=0.0;
480    }
481    solution.setLogLevel(3);
482    solution.dual();
483    double rowObjective[]={1.0,0.5,-10.0};
484    solution.loadProblem(matrix,colLower,colUpper,objective,
485                         rowLower,rowUpper,rowObjective);
486    solution.dual();
487  }
488  {   
489    CoinMpsIO m;
490    std::string fn = mpsDir+"exmip1";
491    m.readMps(fn.c_str(),"mps");
492    ClpSimplex solution;
493    solution.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
494                         m.getObjCoefficients(),
495                         m.getRowLower(),m.getRowUpper());
496    solution.dual();
497    // Test event handling
498    MyEventHandler handler;
499    solution.passInEventHandler(&handler);
500    int numberRows=solution.numberRows();
501    // make sure values pass has something to do
502    for (int i=0;i<numberRows;i++)
503      solution.setRowStatus(i,ClpSimplex::basic);
504    solution.primal(1);
505    assert (solution.secondaryStatus()==102); // Came out at end of pass
506  }
507  // Test Message handler
508  {   
509    CoinMpsIO m;
510    std::string fn = mpsDir+"exmip1";
511    //fn = "Test/subGams4";
512    m.readMps(fn.c_str(),"mps");
513    ClpSimplex model;
514    model.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
515                         m.getObjCoefficients(),
516                         m.getRowLower(),m.getRowUpper());
517    // Message handler
518    MyMessageHandler messageHandler(&model);
519    std::cout<<"Testing derived message handler"<<std::endl;
520    model.passInMessageHandler(&messageHandler);
521    model.messagesPointer()->setDetailMessage(1,102);
522    model.setFactorizationFrequency(10);
523    model.primal();
524
525    // Write saved solutions
526    int nc = model.getNumCols();
527    int s; 
528    std::deque<StdVectorDouble> fep = messageHandler.getFeasibleExtremePoints();
529    int numSavedSolutions = fep.size();
530    for ( s=0; s<numSavedSolutions; ++s ) {
531      const StdVectorDouble & solnVec = fep[s];
532      for ( int c=0; c<nc; ++c ) {
533        if (fabs(solnVec[c])>1.0e-8)
534          std::cout <<"Saved Solution: " <<s <<" ColNum: " <<c <<" Value: " <<solnVec[c] <<std::endl;
535      }
536    }
537    // Solve again without scaling
538    // and maximize then minimize
539    messageHandler.clearFeasibleExtremePoints();
540    model.scaling(0);
541    model.setOptimizationDirection(-1);
542    model.primal();
543    model.setOptimizationDirection(1);
544    model.primal();
545    fep = messageHandler.getFeasibleExtremePoints();
546    numSavedSolutions = fep.size();
547    for ( s=0; s<numSavedSolutions; ++s ) {
548      const StdVectorDouble & solnVec = fep[s];
549      for ( int c=0; c<nc; ++c ) {
550        if (fabs(solnVec[c])>1.0e-8)
551          std::cout <<"Saved Solution: " <<s <<" ColNum: " <<c <<" Value: " <<solnVec[c] <<std::endl;
552      }
553    }
554  }
555  // Test dual ranging
556  {   
557    CoinMpsIO m;
558    std::string fn = mpsDir+"exmip1";
559    m.readMps(fn.c_str(),"mps");
560    ClpSimplex model;
561    model.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
562                         m.getObjCoefficients(),
563                         m.getRowLower(),m.getRowUpper());
564    model.primal();
565    int which[8] = {0,1,2,3,4,5,6,7};
566    double costIncrease[8];
567    int sequenceIncrease[8];
568    double costDecrease[8];
569    int sequenceDecrease[8];
570    // ranging
571    model.dualRanging(8,which,costIncrease,sequenceIncrease,
572                      costDecrease,sequenceDecrease);
573    int i;
574    for ( i=0;i<8;i++)
575      printf("%d increase %g %d, decrease %g %d\n",
576             i,costIncrease[i],sequenceIncrease[i],
577             costDecrease[i],sequenceDecrease[i]);
578    assert (fabs(costDecrease[3])<1.0e-4);
579    assert (fabs(costIncrease[7]-1.0)<1.0e-4);
580    model.setOptimizationDirection(-1);
581    {
582      int j;
583      double * obj = model.objective();
584      int n=model.numberColumns();
585      for (j=0;j<n;j++) 
586        obj[j] *= -1.0;
587    }
588    double costIncrease2[8];
589    int sequenceIncrease2[8];
590    double costDecrease2[8];
591    int sequenceDecrease2[8];
592    // ranging
593    model.dualRanging(8,which,costIncrease2,sequenceIncrease2,
594                      costDecrease2,sequenceDecrease2);
595    for (i=0;i<8;i++) {
596      assert (fabs(costIncrease[i]-costDecrease2[i])<1.0e-6);
597      assert (fabs(costDecrease[i]-costIncrease2[i])<1.0e-6);
598      assert (sequenceIncrease[i]==sequenceDecrease2[i]);
599      assert (sequenceDecrease[i]==sequenceIncrease2[i]);
600    }
601  }
602  // test steepest edge
603  {   
604    CoinMpsIO m;
605    std::string fn = netlibDir+"finnis";
606    m.readMps(fn.c_str(),"mps");
607    ClpModel model;
608    model.loadProblem(*m.getMatrixByCol(),m.getColLower(),
609                    m.getColUpper(),
610                    m.getObjCoefficients(),
611                    m.getRowLower(),m.getRowUpper());
612    ClpSimplex solution(model);
613
614    solution.scaling(1); 
615    solution.setDualBound(1.0e8);
616    //solution.factorization()->maximumPivots(1);
617    //solution.setLogLevel(3);
618    solution.setDualTolerance(1.0e-7);
619    // set objective sense,
620    ClpDualRowSteepest steep;
621    solution.setDualRowPivotAlgorithm(steep);
622    solution.setDblParam(ClpObjOffset,m.objectiveOffset());
623    solution.dual();
624  }
625  // test normal solution
626  {   
627    CoinMpsIO m;
628    std::string fn = netlibDir+"afiro";
629    m.readMps(fn.c_str(),"mps");
630    ClpSimplex solution;
631    ClpModel model;
632    // do twice - without and with scaling
633    int iPass;
634    for (iPass=0;iPass<2;iPass++) {
635      // explicit row objective for testing
636      int nr = m.getNumRows();
637      double * rowObj = new double[nr];
638      CoinFillN(rowObj,nr,0.0);
639      model.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
640                      m.getObjCoefficients(),
641                      m.getRowLower(),m.getRowUpper(),rowObj);
642      delete [] rowObj;
643      solution = ClpSimplex(model);
644      if (iPass) {
645        solution.scaling();
646      }
647      solution.dual();
648      solution.dual();
649      // test optimal
650      assert (solution.status()==0);
651      int numberColumns = solution.numberColumns();
652      int numberRows = solution.numberRows();
653      CoinPackedVector colsol(numberColumns,solution.primalColumnSolution());
654      double * objective = solution.objective();
655      double objValue = colsol.dotProduct(objective);
656      CoinRelFltEq eq(1.0e-8);
657      assert(eq(objValue,-4.6475314286e+02));
658      double * lower = solution.columnLower();
659      double * upper = solution.columnUpper();
660      double * sol = solution.primalColumnSolution();
661      double * result = new double[numberColumns];
662      CoinFillN ( result, numberColumns,0.0);
663      solution.matrix()->transposeTimes(solution.dualRowSolution(), result);
664      int iRow , iColumn;
665      // see if feasible and dual feasible
666      for (iColumn=0;iColumn<numberColumns;iColumn++) {
667        double value = sol[iColumn];
668        assert(value<upper[iColumn]+1.0e-8);
669        assert(value>lower[iColumn]-1.0e-8);
670        value = objective[iColumn]-result[iColumn];
671        assert (value>-1.0e-5);
672        if (sol[iColumn]>1.0e-5)
673          assert (value<1.0e-5);
674      }
675      delete [] result;
676      result = new double[numberRows];
677      CoinFillN ( result, numberRows,0.0);
678      solution.matrix()->times(colsol, result);
679      lower = solution.rowLower();
680      upper = solution.rowUpper();
681      sol = solution.primalRowSolution();
682      for (iRow=0;iRow<numberRows;iRow++) {
683        double value = result[iRow];
684        assert(eq(value,sol[iRow]));
685        assert(value<upper[iRow]+1.0e-8);
686        assert(value>lower[iRow]-1.0e-8);
687      }
688      delete [] result;
689      // test row objective
690      double * rowObjective = solution.rowObjective();
691      CoinDisjointCopyN(solution.dualRowSolution(),numberRows,rowObjective);
692      CoinDisjointCopyN(solution.dualColumnSolution(),numberColumns,objective);
693      // this sets up all slack basis
694      solution.createStatus();
695      solution.dual();
696      CoinFillN(rowObjective,numberRows,0.0);
697      CoinDisjointCopyN(m.getObjCoefficients(),numberColumns,objective);
698      solution.dual();
699    }
700  }
701  // test unbounded
702  {   
703    CoinMpsIO m;
704    std::string fn = netlibDir+"brandy";
705    m.readMps(fn.c_str(),"mps");
706    ClpSimplex solution;
707    // do twice - without and with scaling
708    int iPass;
709    for (iPass=0;iPass<2;iPass++) {
710      solution.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
711                      m.getObjCoefficients(),
712                      m.getRowLower(),m.getRowUpper());
713      if (iPass)
714        solution.scaling();
715      solution.setOptimizationDirection(-1);
716      // test unbounded and ray
717#ifdef DUAL
718      solution.setDualBound(100.0);
719      solution.dual();
720#else
721      solution.primal();
722#endif
723      assert (solution.status()==2);
724      int numberColumns = solution.numberColumns();
725      int numberRows = solution.numberRows();
726      double * lower = solution.columnLower();
727      double * upper = solution.columnUpper();
728      double * sol = solution.primalColumnSolution();
729      double * ray = solution.unboundedRay();
730      double * objective = solution.objective();
731      double objChange=0.0;
732      int iRow , iColumn;
733      // make sure feasible and columns form ray
734      for (iColumn=0;iColumn<numberColumns;iColumn++) {
735        double value = sol[iColumn];
736        assert(value<upper[iColumn]+1.0e-8);
737        assert(value>lower[iColumn]-1.0e-8);
738        value = ray[iColumn];
739        if (value>0.0)
740          assert(upper[iColumn]>1.0e30);
741        else if (value<0.0)
742          assert(lower[iColumn]<-1.0e30);
743        objChange += value*objective[iColumn];
744      }
745      // make sure increasing objective
746      assert(objChange>0.0);
747      double * result = new double[numberRows];
748      CoinFillN ( result, numberRows,0.0);
749      solution.matrix()->times(sol, result);
750      lower = solution.rowLower();
751      upper = solution.rowUpper();
752      sol = solution.primalRowSolution();
753      for (iRow=0;iRow<numberRows;iRow++) {
754        double value = result[iRow];
755        assert(eq(value,sol[iRow]));
756        assert(value<upper[iRow]+2.0e-8);
757        assert(value>lower[iRow]-2.0e-8);
758      }
759      CoinFillN ( result, numberRows,0.0);
760      solution.matrix()->times(ray, result);
761      // there may be small differences (especially if scaled)
762      for (iRow=0;iRow<numberRows;iRow++) {
763        double value = result[iRow];
764        if (value>1.0e-8)
765          assert(upper[iRow]>1.0e30);
766        else if (value<-1.0e-8)
767          assert(lower[iRow]<-1.0e30);
768      }
769      delete [] result;
770      delete [] ray;
771    }
772  }
773  // test infeasible
774  {   
775    CoinMpsIO m;
776    std::string fn = netlibDir+"brandy";
777    m.readMps(fn.c_str(),"mps");
778    ClpSimplex solution;
779    // do twice - without and with scaling
780    int iPass;
781    for (iPass=0;iPass<2;iPass++) {
782      solution.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
783                      m.getObjCoefficients(),
784                      m.getRowLower(),m.getRowUpper());
785      if (iPass)
786        solution.scaling();
787      // test infeasible and ray
788      solution.columnUpper()[0]=0.0;
789#ifdef DUAL
790      solution.setDualBound(100.0);
791      solution.dual();
792#else
793      solution.primal();
794#endif
795      assert (solution.status()==1);
796      int numberColumns = solution.numberColumns();
797      int numberRows = solution.numberRows();
798      double * lower = solution.rowLower();
799      double * upper = solution.rowUpper();
800      double * ray = solution.infeasibilityRay();
801      assert(ray);
802      // construct proof of infeasibility
803      int iRow , iColumn;
804      double lo=0.0,up=0.0;
805      int nl=0,nu=0;
806      for (iRow=0;iRow<numberRows;iRow++) {
807        if (lower[iRow]>-1.0e20) {
808          lo += ray[iRow]*lower[iRow];
809        } else {
810          if (ray[iRow]>1.0e-8) 
811            nl++;
812        }
813        if (upper[iRow]<1.0e20) {
814          up += ray[iRow]*upper[iRow];
815        } else {
816          if (ray[iRow]>1.0e-8) 
817            nu++;
818        }
819      }
820      if (nl)
821        lo=-1.0e100;
822      if (nu)
823        up=1.0e100;
824      double * result = new double[numberColumns];
825      double lo2=0.0,up2=0.0;
826      CoinFillN ( result, numberColumns,0.0);
827      solution.matrix()->transposeTimes(ray, result);
828      lower = solution.columnLower();
829      upper = solution.columnUpper();
830      nl=nu=0;
831      for (iColumn=0;iColumn<numberColumns;iColumn++) {
832        if (result[iColumn]>1.0e-8) {
833          if (lower[iColumn]>-1.0e20)
834            lo2 += result[iColumn]*lower[iColumn];
835          else
836            nl++;
837          if (upper[iColumn]<1.0e20)
838            up2 += result[iColumn]*upper[iColumn];
839          else
840            nu++;
841        } else if (result[iColumn]<-1.0e-8) {
842          if (lower[iColumn]>-1.0e20)
843            up2 += result[iColumn]*lower[iColumn];
844          else
845            nu++;
846          if (upper[iColumn]<1.0e20)
847            lo2 += result[iColumn]*upper[iColumn];
848          else
849            nl++;
850        }
851      }
852      if (nl)
853        lo2=-1.0e100;
854      if (nu)
855        up2=1.0e100;
856      // make sure inconsistency
857      assert(lo2>up||up2<lo);
858      delete [] result;
859      delete [] ray;
860    }
861  }
862  // test delete and add
863  {   
864    CoinMpsIO m;
865    std::string fn = netlibDir+"brandy";
866    m.readMps(fn.c_str(),"mps");
867    ClpSimplex solution;
868    solution.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
869                         m.getObjCoefficients(),
870                         m.getRowLower(),m.getRowUpper());
871    solution.dual();
872    CoinRelFltEq eq(1.0e-8);
873    assert(eq(solution.objectiveValue(),1.5185098965e+03));
874
875    int numberColumns = solution.numberColumns();
876    int numberRows = solution.numberRows();
877    double * saveObj = new double [numberColumns];
878    double * saveLower = new double[numberRows+numberColumns];
879    double * saveUpper = new double[numberRows+numberColumns];
880    int * which = new int [numberRows+numberColumns];
881
882    int numberElements = m.getMatrixByCol()->getNumElements();
883    int * starts = new int[numberRows+numberColumns];
884    int * index = new int[numberElements];
885    double * element = new double[numberElements];
886
887    const CoinBigIndex * startM;
888    const int * lengthM;
889    const int * indexM;
890    const double * elementM;
891
892    int n,nel;
893
894    // delete non basic columns
895    n=0;
896    nel=0;
897    int iRow , iColumn;
898    double * lower = solution.columnLower();
899    double * upper = solution.columnUpper();
900    double * objective = solution.objective();
901    startM = m.getMatrixByCol()->getVectorStarts();
902    lengthM = m.getMatrixByCol()->getVectorLengths();
903    indexM = m.getMatrixByCol()->getIndices();
904    elementM = m.getMatrixByCol()->getElements();
905    starts[0]=0;
906    for (iColumn=0;iColumn<numberColumns;iColumn++) {
907      if (solution.getColumnStatus(iColumn)!=ClpSimplex::basic) {
908        saveObj[n]=objective[iColumn];
909        saveLower[n]=lower[iColumn];
910        saveUpper[n]=upper[iColumn];
911        int j;
912        for (j=startM[iColumn];j<startM[iColumn]+lengthM[iColumn];j++) {
913          index[nel]=indexM[j];
914          element[nel++]=elementM[j];
915        }
916        which[n++]=iColumn;
917        starts[n]=nel;
918      }
919    }
920    solution.deleteColumns(n,which);
921    solution.dual();
922    // Put back
923    solution.addColumns(n,saveLower,saveUpper,saveObj,
924                        starts,index,element);
925    solution.dual();
926    assert(eq(solution.objectiveValue(),1.5185098965e+03));
927
928    // reload with original
929    solution.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
930                         m.getObjCoefficients(),
931                         m.getRowLower(),m.getRowUpper());
932    // delete half rows
933    n=0;
934    nel=0;
935    lower = solution.rowLower();
936    upper = solution.rowUpper();
937    startM = m.getMatrixByRow()->getVectorStarts();
938    lengthM = m.getMatrixByRow()->getVectorLengths();
939    indexM = m.getMatrixByRow()->getIndices();
940    elementM = m.getMatrixByRow()->getElements();
941    starts[0]=0;
942    for (iRow=0;iRow<numberRows;iRow++) {
943      if ((iRow&1)==0) {
944        saveLower[n]=lower[iRow];
945        saveUpper[n]=upper[iRow];
946        int j;
947        for (j=startM[iRow];j<startM[iRow]+lengthM[iRow];j++) {
948          index[nel]=indexM[j];
949          element[nel++]=elementM[j];
950        }
951        which[n++]=iRow;
952        starts[n]=nel;
953      }
954    }
955    solution.deleteRows(n,which);
956    solution.dual();
957    // Put back
958    solution.addRows(n,saveLower,saveUpper,
959                        starts,index,element);
960    solution.dual();
961    assert(eq(solution.objectiveValue(),1.5185098965e+03));
962    // Zero out status array to give some interest
963    memset(solution.statusArray()+numberColumns,0,numberRows);
964    solution.primal(1);
965    assert(eq(solution.objectiveValue(),1.5185098965e+03));
966
967    delete [] saveObj;
968    delete [] saveLower;
969    delete [] saveUpper;
970    delete [] which;
971    delete [] starts;
972    delete [] index;
973    delete [] element;
974  }
975#if 0
976  // Test barrier
977  {
978    CoinMpsIO m;
979    std::string fn = mpsDir+"exmip1";
980    m.readMps(fn.c_str(),"mps");
981    ClpInterior solution;
982    solution.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
983                         m.getObjCoefficients(),
984                         m.getRowLower(),m.getRowUpper());
985    solution.primalDual();
986  }
987#endif
988  // test network
989  //#define QUADRATIC
990#ifndef QUADRATIC
991  if (1) {   
992    std::string fn = mpsDir+"input.130";
993    int numberColumns;
994    int numberRows;
995   
996    FILE * fp = fopen(fn.c_str(),"r");
997    if (!fp) {
998      // Try in Samples
999      fn = "Samples/input.130";
1000      fp = fopen(fn.c_str(),"r");
1001    }
1002    if (!fp) {
1003      fprintf(stderr,"Unable to open file input.130 in mpsDir or Samples directory\n");
1004    } else {
1005      int problem;
1006      char temp[100];
1007      // read and skip
1008      fscanf(fp,"%s",temp);
1009      assert (!strcmp(temp,"BEGIN"));
1010      fscanf(fp,"%*s %*s %d %d %*s %*s %d %*s",&problem, &numberRows, 
1011             &numberColumns);
1012      // scan down to SUPPLY
1013      while (fgets(temp,100,fp)) {
1014        if (!strncmp(temp,"SUPPLY",6))
1015          break;
1016      }
1017      if (strncmp(temp,"SUPPLY",6)) {
1018        fprintf(stderr,"Unable to find SUPPLY\n");
1019        exit(2);
1020      }
1021      // get space for rhs
1022      double * lower = new double[numberRows];
1023      double * upper = new double[numberRows];
1024      int i;
1025      for (i=0;i<numberRows;i++) {
1026        lower[i]=0.0;
1027        upper[i]=0.0;
1028      }
1029      // ***** Remember to convert to C notation
1030      while (fgets(temp,100,fp)) {
1031        int row;
1032        int value;
1033        if (!strncmp(temp,"ARCS",4))
1034          break;
1035        sscanf(temp,"%d %d",&row,&value);
1036        upper[row-1]=-value;
1037        lower[row-1]=-value;
1038      }
1039      if (strncmp(temp,"ARCS",4)) {
1040        fprintf(stderr,"Unable to find ARCS\n");
1041        exit(2);
1042      }
1043      // number of columns may be underestimate
1044      int * head = new int[2*numberColumns];
1045      int * tail = new int[2*numberColumns];
1046      int * ub = new int[2*numberColumns];
1047      int * cost = new int[2*numberColumns];
1048      // ***** Remember to convert to C notation
1049      numberColumns=0;
1050      while (fgets(temp,100,fp)) {
1051        int iHead;
1052        int iTail;
1053        int iUb;
1054        int iCost;
1055        if (!strncmp(temp,"DEMAND",6))
1056          break;
1057        sscanf(temp,"%d %d %d %d",&iHead,&iTail,&iCost,&iUb);
1058        iHead--;
1059        iTail--;
1060        head[numberColumns]=iHead;
1061        tail[numberColumns]=iTail;
1062        ub[numberColumns]=iUb;
1063        cost[numberColumns]=iCost;
1064        numberColumns++;
1065      }
1066      if (strncmp(temp,"DEMAND",6)) {
1067        fprintf(stderr,"Unable to find DEMAND\n");
1068        exit(2);
1069      }
1070      // ***** Remember to convert to C notation
1071      while (fgets(temp,100,fp)) {
1072        int row;
1073        int value;
1074        if (!strncmp(temp,"END",3))
1075          break;
1076        sscanf(temp,"%d %d",&row,&value);
1077        upper[row-1]=value;
1078        lower[row-1]=value;
1079      }
1080      if (strncmp(temp,"END",3)) {
1081        fprintf(stderr,"Unable to find END\n");
1082        exit(2);
1083      }
1084      printf("Problem %d has %d rows and %d columns\n",problem,
1085             numberRows,numberColumns);
1086      fclose(fp);
1087      ClpSimplex  model;
1088      // now build model
1089     
1090      double * objective =new double[numberColumns];
1091      double * lowerColumn = new double[numberColumns];
1092      double * upperColumn = new double[numberColumns];
1093     
1094      double * element = new double [2*numberColumns];
1095      CoinBigIndex * start = new CoinBigIndex [numberColumns+1];
1096      int * row = new int[2*numberColumns];
1097      start[numberColumns]=2*numberColumns;
1098      for (i=0;i<numberColumns;i++) {
1099        start[i]=2*i;
1100        element[2*i]=-1.0;
1101        element[2*i+1]=1.0;
1102        row[2*i]=head[i];
1103        row[2*i+1]=tail[i];
1104        lowerColumn[i]=0.0;
1105        upperColumn[i]=ub[i];
1106        objective[i]=cost[i];
1107      }
1108      // Create Packed Matrix
1109      CoinPackedMatrix matrix;
1110      int * lengths = NULL;
1111      matrix.assignMatrix(true,numberRows,numberColumns,
1112                          2*numberColumns,element,row,start,lengths);
1113      // load model
1114     
1115      model.loadProblem(matrix,
1116                        lowerColumn,upperColumn,objective,
1117                        lower,upper);
1118     
1119      model.factorization()->maximumPivots(200+model.numberRows()/100);
1120      model.createStatus();
1121      double time1 = CoinCpuTime();
1122      model.dual();
1123      std::cout<<"Network problem, ClpPackedMatrix took "<<CoinCpuTime()-time1<<" seconds"<<std::endl;
1124      ClpPlusMinusOneMatrix plusMinus(matrix);
1125      assert (plusMinus.getIndices()); // would be zero if not +- one
1126      model.loadProblem(plusMinus,
1127                        lowerColumn,upperColumn,objective,
1128                        lower,upper);
1129     
1130      model.factorization()->maximumPivots(200+model.numberRows()/100);
1131      model.createStatus();
1132      time1 = CoinCpuTime();
1133      model.dual();
1134      std::cout<<"Network problem, ClpPlusMinusOneMatrix took "<<CoinCpuTime()-time1<<" seconds"<<std::endl;
1135      ClpNetworkMatrix network(numberColumns,head,tail);
1136      model.loadProblem(network,
1137                        lowerColumn,upperColumn,objective,
1138                        lower,upper);
1139     
1140      model.factorization()->maximumPivots(200+model.numberRows()/100);
1141      model.createStatus();
1142      time1 = CoinCpuTime();
1143      model.dual();
1144      std::cout<<"Network problem, ClpNetworkMatrix took "<<CoinCpuTime()-time1<<" seconds"<<std::endl;
1145      delete [] lower;
1146      delete [] upper;
1147      delete [] head;
1148      delete [] tail;
1149      delete [] ub;
1150      delete [] cost;
1151      delete [] objective;
1152      delete [] lowerColumn;
1153      delete [] upperColumn;
1154    }
1155  }
1156#endif
1157#ifdef QUADRATIC
1158  // Test quadratic to solve linear
1159  if (1) {   
1160    CoinMpsIO m;
1161    std::string fn = mpsDir+"exmip1";
1162    m.readMps(fn.c_str(),"mps");
1163    ClpSimplex solution;
1164    solution.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
1165                         m.getObjCoefficients(),
1166                         m.getRowLower(),m.getRowUpper());
1167    //solution.dual();
1168    // get quadratic part
1169    int numberColumns=solution.numberColumns();
1170    int * start=new int [numberColumns+1];
1171    int * column = new int[numberColumns];
1172    double * element = new double[numberColumns];
1173    int i;
1174    start[0]=0;
1175    int n=0;
1176    int kk=numberColumns-1;
1177    int kk2=numberColumns-1;
1178    for (i=0;i<numberColumns;i++) {
1179      if (i>=kk) {
1180        column[n]=i;
1181        if (i>=kk2)
1182          element[n]=1.0e-1;
1183        else
1184          element[n]=0.0;
1185        n++;
1186      }
1187      start[i+1]=n;
1188    }
1189    // Load up objective
1190    solution.loadQuadraticObjective(numberColumns,start,column,element);
1191    delete [] start;
1192    delete [] column;
1193    delete [] element;
1194    //solution.quadraticSLP(50,1.0e-4);
1195    double objValue = solution.getObjValue();
1196    CoinRelFltEq eq(1.0e-4);
1197    //assert(eq(objValue,820.0));
1198    solution.setLogLevel(63);
1199    solution.quadraticPrimal(1);
1200    objValue = solution.getObjValue();
1201    //assert(eq(objValue,3.2368421));
1202    //exit(77);
1203  }
1204  // Test quadratic
1205  if (1) {   
1206    CoinMpsIO m;
1207    std::string fn = mpsDir+"share2qp";
1208    //fn = "share2qpb";
1209    m.readMps(fn.c_str(),"mps");
1210    ClpSimplex model;
1211    model.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
1212                         m.getObjCoefficients(),
1213                         m.getRowLower(),m.getRowUpper());
1214    model.dual();
1215    // get quadratic part
1216    int * start=NULL;
1217    int * column = NULL;
1218    double * element = NULL;
1219    m.readQuadraticMps(NULL,start,column,element,2);
1220    int column2[200];
1221    double element2[200];
1222    int start2[80];
1223    int j;
1224    start2[0]=0;
1225    int nel=0;
1226    bool good=false;
1227    for (j=0;j<79;j++) {
1228      if (start[j]==start[j+1]) {
1229        column2[nel]=j;
1230        element2[nel]=0.0;
1231        nel++;
1232      } else {
1233        int i;
1234        for (i=start[j];i<start[j+1];i++) {
1235          column2[nel]=column[i];
1236          element2[nel++]=element[i];
1237        }
1238      }
1239      start2[j+1]=nel;
1240    }
1241    // Load up objective
1242    if (good)
1243      model.loadQuadraticObjective(model.numberColumns(),start2,column2,element2);
1244    else
1245      model.loadQuadraticObjective(model.numberColumns(),start,column,element);
1246    delete [] start;
1247    delete [] column;
1248    delete [] element;
1249    int numberColumns=model.numberColumns();
1250#if 0
1251    model.quadraticSLP(50,1.0e-4);
1252#else
1253    // Get feasible
1254    ClpObjective * saveObjective = model.objectiveAsObject()->clone();
1255    ClpLinearObjective zeroObjective(NULL,numberColumns);
1256    model.setObjective(&zeroObjective);
1257    model.dual();
1258    model.setObjective(saveObjective);
1259    delete saveObjective;
1260#endif
1261    model.setLogLevel(63);
1262    //exit(77);
1263    model.setFactorizationFrequency(10);
1264    model.quadraticPrimal(1);
1265    double objValue = model.getObjValue();
1266    const double * solution = model.primalColumnSolution();
1267    int i;
1268    for (i=0;i<numberColumns;i++)
1269      if (solution[i])
1270        printf("%d %g\n",i,solution[i]);
1271    CoinRelFltEq eq(1.0e-4);
1272    assert(eq(objValue,-400.92));
1273  }
1274  if (0) {   
1275    CoinMpsIO m;
1276    std::string fn = "./beale";
1277    //fn = "./jensen";
1278    m.readMps(fn.c_str(),"mps");
1279    ClpSimplex solution;
1280    solution.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
1281                         m.getObjCoefficients(),
1282                         m.getRowLower(),m.getRowUpper());
1283    solution.setDblParam(ClpObjOffset,m.objectiveOffset());
1284    solution.dual();
1285    // get quadratic part
1286    int * start=NULL;
1287    int * column = NULL;
1288    double * element = NULL;
1289    m.readQuadraticMps(NULL,start,column,element,2);
1290    // Load up objective
1291    solution.loadQuadraticObjective(solution.numberColumns(),start,column,element);
1292    delete [] start;
1293    delete [] column;
1294    delete [] element;
1295    solution.quadraticPrimal(1);
1296    solution.quadraticSLP(50,1.0e-4);
1297    double objValue = solution.getObjValue();
1298    CoinRelFltEq eq(1.0e-4);
1299    assert(eq(objValue,0.5));
1300    solution.quadraticPrimal();
1301    objValue = solution.getObjValue();
1302    assert(eq(objValue,0.5));
1303  }
1304#endif 
1305}
Note: See TracBrowser for help on using the repository browser.