source: branches/devel-1/Test/unitTest.cpp @ 29

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

Presolve (no changes to Makefile)

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 33.3 KB
Line 
1// Copyright (C) 2002, International Business Machines
2// Corporation and others.  All Rights Reserved.
3#if defined(_MSC_VER)
4// Turn off compiler warning about long names
5#  pragma warning(disable:4786)
6#endif
7
8#include <cassert>
9#include <cstdio>
10#include <cmath>
11#include <cfloat>
12#include <string>
13#include <iostream>
14
15#include <time.h>
16#include <sys/times.h>
17#include <sys/resource.h>
18#include <unistd.h>
19
20#include "CoinMpsIO.hpp"
21#include "CoinPackedMatrix.hpp"
22#include "CoinPackedVector.hpp"
23#include "CoinHelperFunctions.hpp"
24
25#include "ClpFactorization.hpp"
26#include "ClpSimplex.hpp"
27#include "ClpDualRowSteepest.hpp"
28#include "ClpDualRowDantzig.hpp"
29#include "ClpPrimalColumnSteepest.hpp"
30#include "ClpPrimalColumnDantzig.hpp"
31#include "ClpParameters.hpp"
32
33#include "Presolve.hpp"
34
35//#############################################################################
36
37#ifdef NDEBUG
38#undef NDEBUG
39#endif
40
41// Function Prototypes. Function definitions is in this file.
42void testingMessage( const char * const msg );
43
44//----------------------------------------------------------------
45// unitTest [-mpsDir=V1] [-netlibDir=V2] [-netlib]
46//
47// where:
48//   -mpsDir: directory containing mps test files
49//       Default value V1="../Mps/Sample"   
50//   -netlibDir: directory containing netlib files
51//       Default value V2="../Mps/Netlib"
52//   -netlib
53//       If specified, then netlib test set run
54//
55// All parameters are optional.
56//----------------------------------------------------------------
57
58int mainTest (int argc, const char *argv[],bool doDual,
59              ClpSimplex empty, bool doPresolve)
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
144#if 1
145    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);
146    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);
147    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);
148   
149    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);
150    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);
151    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);
152    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);
153    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);
154    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);
155    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);
156    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);
157    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);
158    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);
159    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);
160    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);
161    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);
162    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);
163    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);
164    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);
165    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);
166    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);
167    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);
168    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);
169    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);
170    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);
171    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);
172    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);
173    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);
174    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.
175    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 );
176    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);
177    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);
178    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);
179    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);
180    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);
181    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);
182    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);
183    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);
184    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);
185#endif
186    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);
187    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);
188    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);
189    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);
190    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);
191    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);
192    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);
193    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);
194    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);
195    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);
196    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);
197    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);
198    //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);
199    //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);
200    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);
201    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);
202    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);
203    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);
204    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);
205    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);
206    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);
207    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);
208    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);
209    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);
210    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);
211    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);
212    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);
213    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);
214    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);
215    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);
216    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);
217    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);
218    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);
219    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);
220    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);
221    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);
222    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);
223    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);
224    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);
225    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);
226    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);
227    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);
228    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);
229    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);
230    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);
231    //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);
232    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); 
233    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);
234    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);
235    //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);
236    //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);
237    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);
238    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);
239    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);
240    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);
241  // Loop once for each Mps File
242    for (m=0; m<mpsName.size(); m++ ) {
243      std::cerr <<"  processing mps file: " <<mpsName[m] 
244                <<" (" <<m+1 <<" out of " <<mpsName.size() <<")" <<std::endl;
245   
246      // Read data mps file,
247      std::string fn = netlibDir+mpsName[m];
248      CoinMpsIO mps;
249      mps.readMps(fn.c_str(),"mps");
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 (doPresolve) {
258        Presolve pinfo;
259        ClpSimplex * model2 = pinfo.presolvedModel(solution,1.0e-8);
260        if (doDual)
261          model2->dual();
262        else
263          model2->primal();
264        pinfo.postsolve(true);
265       
266        delete model2;
267        printf("Resolving from postsolved model\n");
268        // later try without (1) and check duals before solve
269        solution.primal(1);
270        if (solution.numberIterations())
271          printf("****** iterated %d\n",solution.numberIterations());
272        solution.checkSolution();
273        printf("%g dual %g(%d) Primal %g(%d)\n",
274               solution.objectiveValue(),
275               solution.sumDualInfeasibilities(),
276               solution.numberDualInfeasibilities(),
277               solution.sumPrimalInfeasibilities(),
278               solution.numberPrimalInfeasibilities());
279        {
280          Presolve pinfoA;
281          model2 = pinfoA.presolvedModel(solution,1.0e-8);
282
283          printf("Resolving from presolved optimal solution\n");
284          model2->primal(1);
285               
286          delete model2;
287        }
288
289      } else {
290        if (doDual) {
291          solution.dual();
292        } else {
293          solution.primal();
294        }
295      }
296      // Test objective solution value
297      {
298        double soln = solution.objectiveValue();
299        CoinRelFltEq eq(objValueTol[m]);
300        std::cerr <<soln <<",  " <<objValue[m] <<" diff "<<
301          soln-objValue[m]<<std::endl;
302        if(!eq(soln,objValue[m]))
303          printf("** difference fails\n");
304      }
305    }
306   
307  }
308  else {
309    testingMessage( "***Skipped Testing on netlib    ***\n" );
310    testingMessage( "***use -netlib to test class***\n" );
311  }
312
313  testingMessage( "All tests completed successfully\n" );
314  return 0;
315}
316
317 
318// Display message on stdout and stderr
319void testingMessage( const char * const msg )
320{
321  std::cerr <<msg;
322  //cout <<endl <<"*****************************************"
323  //     <<endl <<msg <<endl;
324}
325
326//--------------------------------------------------------------------------
327// test factorization methods and simplex method
328void
329ClpSimplexUnitTest(const std::string & mpsDir,
330                   const std::string & netlibDir)
331{
332 
333  CoinRelFltEq eq(0.000001);
334
335  {
336    ClpSimplex solution;
337 
338    // matrix data
339    //deliberate hiccup of 2 between 0 and 1
340    int start[5]={0,4,7,8,9};
341    int length[5]={2,3,1,1,1};
342    int rows[11]={0,2,-1,-1,0,1,2,0,1,2};
343    double elements[11]={7.0,2.0,1.0e10,1.0e10,-2.0,1.0,-2.0,1,1,1};
344    CoinPackedMatrix matrix(true,3,5,8,elements,rows,start,length);
345   
346    // rim data
347    double objective[7]={-4.0,1.0,0.0,0.0,0.0,0.0,0.0};
348    double rowLower[5]={14.0,3.0,3.0,1.0e10,1.0e10};
349    double rowUpper[5]={14.0,3.0,3.0,-1.0e10,-1.0e10};
350    double colLower[7]={0.0,0.0,0.0,0.0,0.0,0.0,0.0};
351    double colUpper[7]={100.0,100.0,100.0,100.0,100.0,100.0,100.0};
352   
353    // basis 1
354    int rowBasis1[3]={-1,-1,-1};
355    int colBasis1[5]={1,1,-1,-1,1};
356    solution.loadProblem(matrix,colLower,colUpper,objective,
357                         rowLower,rowUpper);
358    int i;
359    solution.createStatus();
360    for (i=0;i<3;i++) {
361      if (rowBasis1[i]<0) {
362        solution.setRowStatus(i,ClpSimplex::atLowerBound);
363      } else {
364        solution.setRowStatus(i,ClpSimplex::basic);
365      }
366    }
367    for (i=0;i<5;i++) {
368      if (colBasis1[i]<0) {
369        solution.setColumnStatus(i,ClpSimplex::atLowerBound);
370      } else {
371        solution.setColumnStatus(i,ClpSimplex::basic);
372      }
373    }
374    solution.setLogLevel(3+4+8+16+32);
375    solution.primal();
376    for (i=0;i<3;i++) {
377      if (rowBasis1[i]<0) {
378        solution.setRowStatus(i,ClpSimplex::atLowerBound);
379      } else {
380        solution.setRowStatus(i,ClpSimplex::basic);
381      }
382    }
383    for (i=0;i<5;i++) {
384      if (colBasis1[i]<0) {
385        solution.setColumnStatus(i,ClpSimplex::atLowerBound);
386      } else {
387        solution.setColumnStatus(i,ClpSimplex::basic);
388      }
389    }
390    // intricate stuff does not work with scaling
391    solution.scaling(0);
392    assert(!solution.factorize ( ));
393    const double * colsol = solution.primalColumnSolution();
394    const double * rowsol = solution.primalRowSolution();
395    solution.getSolution(rowsol,colsol);
396    double colsol1[5]={20.0/7.0,3.0,0.0,0.0,23.0/7.0};
397    for (i=0;i<5;i++) {
398      assert(eq(colsol[i],colsol1[i]));
399    }
400    // now feed in again without actually doing factorization
401    ClpFactorization factorization2 = *solution.factorization();
402    ClpSimplex solution2 = solution;
403    solution2.setFactorization(factorization2);
404    solution2.createStatus();
405    for (i=0;i<3;i++) {
406      if (rowBasis1[i]<0) {
407        solution2.setRowStatus(i,ClpSimplex::atLowerBound);
408      } else {
409        solution2.setRowStatus(i,ClpSimplex::basic);
410      }
411    }
412    for (i=0;i<5;i++) {
413      if (colBasis1[i]<0) {
414        solution2.setColumnStatus(i,ClpSimplex::atLowerBound);
415      } else {
416        solution2.setColumnStatus(i,ClpSimplex::basic);
417      }
418    }
419    // intricate stuff does not work with scaling
420    solution2.scaling(0);
421    solution2.getSolution(rowsol,colsol);
422    colsol = solution2.primalColumnSolution();
423    rowsol = solution2.primalRowSolution();
424    for (i=0;i<5;i++) {
425      assert(eq(colsol[i],colsol1[i]));
426    }
427    solution2.setDualBound(0.1);
428    solution2.dual();
429    objective[2]=-1.0;
430    objective[3]=-0.5;
431    objective[4]=10.0;
432    solution.dual();
433    for (i=0;i<3;i++) {
434      rowLower[i]=-1.0e50;
435      colUpper[i+2]=0.0;
436    }
437    solution.setLogLevel(3);
438    solution.dual();
439    double rowObjective[]={1.0,0.5,-10.0};
440    solution.loadProblem(matrix,colLower,colUpper,objective,
441                         rowLower,rowUpper,rowObjective);
442    solution.dual();
443  }
444  {   
445    CoinMpsIO m;
446    std::string fn = mpsDir+"exmip1";
447    m.readMps(fn.c_str(),"mps");
448    ClpSimplex solution;
449    solution.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
450                         m.getObjCoefficients(),
451                         m.getRowLower(),m.getRowUpper());
452    solution.dual();
453  }
454  // test steepest edge
455  {   
456    CoinMpsIO m;
457    std::string fn = netlibDir+"finnis";
458    m.readMps(fn.c_str(),"mps");
459    ClpModel model;
460    model.loadProblem(*m.getMatrixByCol(),m.getColLower(),
461                    m.getColUpper(),
462                    m.getObjCoefficients(),
463                    m.getRowLower(),m.getRowUpper());
464    ClpSimplex solution(model);
465
466    solution.scaling(1); 
467    solution.setDualBound(1.0e8);
468    //solution.factorization()->maximumPivots(1);
469    //solution.setLogLevel(3);
470    solution.setDualTolerance(1.0e-7);
471    // set objective sense,
472    ClpDualRowSteepest steep;
473    solution.setDualRowPivotAlgorithm(steep);
474    solution.setDblParam(ClpObjOffset,m.objectiveOffset());
475    solution.dual();
476  }
477  // test normal solution
478  {   
479    CoinMpsIO m;
480    std::string fn = netlibDir+"afiro";
481    m.readMps(fn.c_str(),"mps");
482    ClpSimplex solution;
483    ClpModel model;
484    // do twice - without and with scaling
485    int iPass;
486    for (iPass=0;iPass<2;iPass++) {
487      // explicit row objective for testing
488      int nr = m.getNumRows();
489      double * rowObj = new double[nr];
490      CoinFillN(rowObj,nr,0.0);
491      model.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
492                      m.getObjCoefficients(),
493                      m.getRowLower(),m.getRowUpper(),rowObj);
494      delete [] rowObj;
495      solution = ClpSimplex(model);
496      if (iPass) {
497        solution.scaling();
498      }
499      solution.dual();
500      solution.dual();
501      // test optimal
502      assert (solution.status()==0);
503      int numberColumns = solution.numberColumns();
504      int numberRows = solution.numberRows();
505      CoinPackedVector colsol(numberColumns,solution.primalColumnSolution());
506      double * objective = solution.objective();
507      double objValue = colsol.dotProduct(objective);
508      CoinRelFltEq eq(1.0e-8);
509      assert(eq(objValue,-4.6475314286e+02));
510      double * lower = solution.columnLower();
511      double * upper = solution.columnUpper();
512      double * sol = solution.primalColumnSolution();
513      double * result = new double[numberColumns];
514      CoinFillN ( result, numberColumns,0.0);
515      solution.matrix()->transposeTimes(solution.dualRowSolution(), result);
516      int iRow , iColumn;
517      // see if feasible and dual feasible
518      for (iColumn=0;iColumn<numberColumns;iColumn++) {
519        double value = sol[iColumn];
520        assert(value<upper[iColumn]+1.0e-8);
521        assert(value>lower[iColumn]-1.0e-8);
522        value = objective[iColumn]-result[iColumn];
523        assert (value>-1.0e-5);
524        if (sol[iColumn]>1.0e-5)
525          assert (value<1.0e-5);
526      }
527      delete [] result;
528      result = new double[numberRows];
529      CoinFillN ( result, numberRows,0.0);
530      solution.matrix()->times(colsol, result);
531      lower = solution.rowLower();
532      upper = solution.rowUpper();
533      sol = solution.primalRowSolution();
534      for (iRow=0;iRow<numberRows;iRow++) {
535        double value = result[iRow];
536        assert(eq(value,sol[iRow]));
537        assert(value<upper[iRow]+1.0e-8);
538        assert(value>lower[iRow]-1.0e-8);
539      }
540      delete [] result;
541      // test row objective
542      double * rowObjective = solution.rowObjective();
543      CoinDisjointCopyN(solution.dualRowSolution(),numberRows,rowObjective);
544      CoinDisjointCopyN(solution.dualColumnSolution(),numberColumns,objective);
545      // this sets up all slack basis
546      solution.createStatus();
547      solution.dual();
548      CoinFillN(rowObjective,numberRows,0.0);
549      CoinDisjointCopyN(m.getObjCoefficients(),numberColumns,objective);
550      solution.dual();
551    }
552  }
553  // test unbounded
554  {   
555    CoinMpsIO m;
556    std::string fn = netlibDir+"brandy";
557    m.readMps(fn.c_str(),"mps");
558    ClpSimplex solution;
559    // do twice - without and with scaling
560    int iPass;
561    for (iPass=0;iPass<2;iPass++) {
562      solution.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
563                      m.getObjCoefficients(),
564                      m.getRowLower(),m.getRowUpper());
565      if (iPass)
566        solution.scaling();
567      solution.setOptimizationDirection(-1);
568      // test unbounded and ray
569#ifdef DUAL
570      solution.setDualBound(100.0);
571      solution.dual();
572#else
573      solution.primal();
574#endif
575      assert (solution.status()==2);
576      int numberColumns = solution.numberColumns();
577      int numberRows = solution.numberRows();
578      double * lower = solution.columnLower();
579      double * upper = solution.columnUpper();
580      double * sol = solution.primalColumnSolution();
581      double * ray = solution.unboundedRay();
582      double * objective = solution.objective();
583      double objChange=0.0;
584      int iRow , iColumn;
585      // make sure feasible and columns form ray
586      for (iColumn=0;iColumn<numberColumns;iColumn++) {
587        double value = sol[iColumn];
588        assert(value<upper[iColumn]+1.0e-8);
589        assert(value>lower[iColumn]-1.0e-8);
590        value = ray[iColumn];
591        if (value>0.0)
592          assert(upper[iColumn]>1.0e30);
593        else if (value<0.0)
594          assert(lower[iColumn]<-1.0e30);
595        objChange += value*objective[iColumn];
596      }
597      // make sure increasing objective
598      assert(objChange>0.0);
599      double * result = new double[numberRows];
600      CoinFillN ( result, numberRows,0.0);
601      solution.matrix()->times(sol, result);
602      lower = solution.rowLower();
603      upper = solution.rowUpper();
604      sol = solution.primalRowSolution();
605      for (iRow=0;iRow<numberRows;iRow++) {
606        double value = result[iRow];
607        assert(eq(value,sol[iRow]));
608        assert(value<upper[iRow]+1.0e-8);
609        assert(value>lower[iRow]-1.0e-8);
610      }
611      CoinFillN ( result, numberRows,0.0);
612      solution.matrix()->times(ray, result);
613      // there may be small differences (especially if scaled)
614      for (iRow=0;iRow<numberRows;iRow++) {
615        double value = result[iRow];
616        if (value>1.0e-8)
617          assert(upper[iRow]>1.0e30);
618        else if (value<-1.0e-8)
619          assert(lower[iRow]<-1.0e30);
620      }
621      delete [] result;
622      delete [] ray;
623    }
624  }
625  // test infeasible
626  {   
627    CoinMpsIO m;
628    std::string fn = netlibDir+"brandy";
629    m.readMps(fn.c_str(),"mps");
630    ClpSimplex solution;
631    // do twice - without and with scaling
632    int iPass;
633    for (iPass=0;iPass<2;iPass++) {
634      solution.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
635                      m.getObjCoefficients(),
636                      m.getRowLower(),m.getRowUpper());
637      if (iPass)
638        solution.scaling();
639      // test infeasible and ray
640      solution.columnUpper()[0]=0.0;
641#ifdef DUAL
642      solution.setDualBound(100.0);
643      solution.dual();
644#else
645      solution.primal();
646#endif
647      assert (solution.status()==1);
648      int numberColumns = solution.numberColumns();
649      int numberRows = solution.numberRows();
650      double * lower = solution.rowLower();
651      double * upper = solution.rowUpper();
652      double * ray = solution.infeasibilityRay();
653      assert(ray);
654      // construct proof of infeasibility
655      int iRow , iColumn;
656      double lo=0.0,up=0.0;
657      int nl=0,nu=0;
658      for (iRow=0;iRow<numberRows;iRow++) {
659        if (lower[iRow]>-1.0e20) {
660          lo += ray[iRow]*lower[iRow];
661        } else {
662          if (ray[iRow]>1.0e-8) 
663            nl++;
664        }
665        if (upper[iRow]<1.0e20) {
666          up += ray[iRow]*upper[iRow];
667        } else {
668          if (ray[iRow]>1.0e-8) 
669            nu++;
670        }
671      }
672      if (nl)
673        lo=-1.0e100;
674      if (nu)
675        up=1.0e100;
676      double * result = new double[numberColumns];
677      double lo2=0.0,up2=0.0;
678      CoinFillN ( result, numberColumns,0.0);
679      solution.matrix()->transposeTimes(ray, result);
680      lower = solution.columnLower();
681      upper = solution.columnUpper();
682      nl=nu=0;
683      for (iColumn=0;iColumn<numberColumns;iColumn++) {
684        if (result[iColumn]>1.0e-8) {
685          if (lower[iColumn]>-1.0e20)
686            lo2 += result[iColumn]*lower[iColumn];
687          else
688            nl++;
689          if (upper[iColumn]<1.0e20)
690            up2 += result[iColumn]*upper[iColumn];
691          else
692            nu++;
693        } else if (result[iColumn]<-1.0e-8) {
694          if (lower[iColumn]>-1.0e20)
695            up2 += result[iColumn]*lower[iColumn];
696          else
697            nu++;
698          if (upper[iColumn]<1.0e20)
699            lo2 += result[iColumn]*upper[iColumn];
700          else
701            nl++;
702        }
703      }
704      if (nl)
705        lo2=-1.0e100;
706      if (nu)
707        up2=1.0e100;
708      // make sure inconsistency
709      assert(lo2>up||up2<lo);
710      delete [] result;
711      delete [] ray;
712    }
713  }
714 
715}
Note: See TracBrowser for help on using the repository browser.