CLP: unitTest.cpp Source File

00001 // Copyright (C) 2002, International Business Machines
00002 // Corporation and others.  All Rights Reserved.
00003 
00004 #include "CoinPragma.hpp"
00005 #include <cassert>
00006 #include <cstdio>
00007 #include <cmath>
00008 #include <cfloat>
00009 #include <string>
00010 #include <iostream>
00011 
00012 #include "CoinMpsIO.hpp"
00013 #include "CoinPackedMatrix.hpp"
00014 #include "CoinPackedVector.hpp"
00015 #include "CoinHelperFunctions.hpp"
00016 #include "CoinTime.hpp"
00017 
00018 #include "ClpFactorization.hpp"
00019 #include "ClpSimplex.hpp"
00020 //#include "ClpInterior.hpp"
00021 #include "ClpLinearObjective.hpp"
00022 #include "ClpDualRowSteepest.hpp"
00023 #include "ClpDualRowDantzig.hpp"
00024 #include "ClpPrimalColumnSteepest.hpp"
00025 #include "ClpPrimalColumnDantzig.hpp"
00026 #include "ClpParameters.hpp"
00027 #include "ClpNetworkMatrix.hpp"
00028 #include "ClpPlusMinusOneMatrix.hpp"
00029 #include "MyMessageHandler.hpp"
00030 
00031 #include "ClpPresolve.hpp"
00032 #include "Idiot.hpp"
00033 
00034 
00035 //#############################################################################
00036 
00037 #ifdef NDEBUG
00038 #undef NDEBUG
00039 #endif
00040 
00041 // Function Prototypes. Function definitions is in this file.
00042 void testingMessage( const char * const msg );
00043 
00044 //----------------------------------------------------------------
00045 // unitTest [-mpsDir=V1] [-netlibDir=V2] [-netlib]
00046 // 
00047 // where:
00048 //   -mpsDir: directory containing mps test files
00049 //       Default value V1="../Mps/Sample"    
00050 //   -netlibDir: directory containing netlib files
00051 //       Default value V2="../Mps/Netlib"
00052 //   -netlib
00053 //       If specified, then netlib test set run
00054 //
00055 // All parameters are optional.
00056 //----------------------------------------------------------------
00057 int mainTest (int argc, const char *argv[],bool doDual,
00058               ClpSimplex empty, bool doPresolve,int doIdiot)
00059 {
00060   int i;
00061 
00062   // define valid parameter keywords
00063   std::set<std::string> definedKeyWords;
00064   definedKeyWords.insert("-mpsDir");
00065   definedKeyWords.insert("-netlibDir");
00066   definedKeyWords.insert("-netlib");
00067 
00068   // Create a map of parameter keys and associated data
00069   std::map<std::string,std::string> parms;
00070   for ( i=1; i<argc; i++ ) {
00071     std::string parm(argv[i]);
00072     std::string key,value;
00073     unsigned int  eqPos = parm.find('=');
00074 
00075     // Does parm contain and '='
00076     if ( eqPos==std::string::npos ) {
00077       //Parm does not contain '='
00078       key = parm;
00079     }
00080     else {
00081       key=parm.substr(0,eqPos);
00082       value=parm.substr(eqPos+1);
00083     }
00084 
00085     // Is specifed key valid?
00086     if ( definedKeyWords.find(key) == definedKeyWords.end() ) {
00087       // invalid key word.
00088       // Write help text
00089       std::cerr <<"Undefined parameter \"" <<key <<"\".\n";
00090       std::cerr <<"Correct usage: \n";
00091       std::cerr <<"  unitTest [-mpsDir=V1] [-netlibDir=V2] [-netlib]\n";
00092       std::cerr <<"  where:\n";
00093       std::cerr <<"    -mpsDir: directory containing mps test files\n";
00094       std::cerr <<"        Default value V1=\"../Mps/Sample\"\n";
00095       std::cerr <<"    -netlibDir: directory containing netlib files\n";
00096       std::cerr <<"        Default value V2=\"../Mps/Netlib\"\n";
00097       std::cerr <<"    -netlib\n";
00098       std::cerr <<"        If specified, then netlib testset run.\n";
00099       return 1;
00100     }
00101     parms[key]=value;
00102   }
00103   
00104   const char dirsep =  CoinFindDirSeparator();
00105   // Set directory containing mps data files.
00106   std::string mpsDir;
00107   if (parms.find("-mpsDir") != parms.end())
00108     mpsDir=parms["-mpsDir"] + dirsep;
00109   else 
00110     mpsDir = dirsep == '/' ? "../Mps/Sample/" : "..\\Mps\\Sample\\";
00111  
00112   // Set directory containing netlib data files.
00113   std::string netlibDir;
00114   if (parms.find("-netlibDir") != parms.end())
00115     netlibDir=parms["-netlibDir"] + dirsep;
00116   else 
00117     netlibDir = dirsep == '/' ? "../Mps/Netlib/" : "..\\Mps\\Netlib\\";
00118 
00119   testingMessage( "Testing ClpSimplex\n" );
00120   ClpSimplexUnitTest(mpsDir,netlibDir);
00121   if (parms.find("-netlib") != parms.end())
00122   {
00123     unsigned int m;
00124     
00125     // Define test problems: 
00126     //   mps names, 
00127     //   maximization or minimization, 
00128     //   Number of rows and columns in problem, and
00129     //   objective function value
00130     std::vector<std::string> mpsName;
00131     std::vector<bool> min;
00132     std::vector<int> nRows;
00133     std::vector<int> nCols;
00134     std::vector<double> objValue;
00135     std::vector<double> objValueTol;
00136     mpsName.push_back("25fv47");
00137     min.push_back(true);
00138     nRows.push_back(822);
00139     nCols.push_back(1571);
00140     objValueTol.push_back(1.E-10);
00141     objValue.push_back(5.5018458883E+03);
00142     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);
00143     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);
00144     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);
00145     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);
00146     
00147     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);
00148     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);
00149     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);
00150     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);
00151     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);
00152     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);
00153     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);
00154     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);
00155     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);
00156     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);
00157     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);
00158     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);
00159     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);
00160     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);
00161     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);
00162     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);
00163     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);
00164     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);
00165     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);
00166     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);
00167     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);
00168     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);
00169     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);
00170     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);
00171     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.
00172     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 );
00173     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);
00174     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);
00175     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);
00176     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);
00177     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);
00178     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);
00179     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);
00180     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);
00181     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);
00182     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);
00183     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);
00184     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);
00185     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);
00186     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);
00187     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);
00188     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);
00189     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);
00190     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);
00191     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);
00192     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);
00193     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);
00194     //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);
00195     //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);
00196     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);
00197     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);
00198     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);
00199     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);
00200     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);
00201     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);
00202     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);
00203     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);
00204     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);
00205     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);
00206     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);
00207     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);
00208     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);
00209     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);
00210     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);
00211     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);
00212     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);
00213     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);
00214     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);
00215     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);
00216     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);
00217     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);
00218     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);
00219     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);
00220     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);
00221     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);
00222     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);
00223     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);
00224     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);
00225     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);
00226     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);
00227     //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); 
00228     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); 
00229     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);
00230     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);
00231     //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);
00232     //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);
00233     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);
00234     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);
00235     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);
00236     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);
00237 
00238     double timeTaken =0.0;
00239   // Loop once for each Mps File
00240     for (m=0; m<mpsName.size(); m++ ) {
00241       std::cerr <<"  processing mps file: " <<mpsName[m] 
00242                 <<" (" <<m+1 <<" out of " <<mpsName.size() <<")" <<std::endl;
00243     
00244       // Read data mps file,
00245       std::string fn = netlibDir+mpsName[m];
00246       CoinMpsIO mps;
00247       mps.readMps(fn.c_str(),"mps");
00248       double time1 = CoinCpuTime();
00249       ClpSimplex solution=empty;
00250       solution.loadProblem(*mps.getMatrixByCol(),mps.getColLower(),
00251                            mps.getColUpper(),
00252                            mps.getObjCoefficients(),
00253                            mps.getRowLower(),mps.getRowUpper());
00254 
00255       solution.setDblParam(ClpObjOffset,mps.objectiveOffset());
00256 #if 0
00257       solution.setOptimizationDirection(-1);
00258       {
00259         int j;
00260         double * obj = solution.objective();
00261         int n=solution.numberColumns();
00262         for (j=0;j<n;j++) 
00263           obj[j] *= -1.0;
00264       }
00265 #endif
00266       if (doPresolve) {
00267         ClpPresolve pinfo;
00268         double a=CoinCpuTime();
00269         ClpSimplex * model2 = pinfo.presolvedModel(solution,1.0e-8,
00270                                                    false,5,true);
00271         printf("Presolve took %g seconds\n",CoinCpuTime()-a);
00272         // change from 200 (unless user has changed)
00273         if (model2->factorization()->maximumPivots()==200&&
00274             model2->numberRows()>-5000)
00275           model2->factorization()->maximumPivots(100+model2->numberRows()/100);
00276         if (doDual) {
00277           // faster if bounds tightened
00278           int numberInfeasibilities = model2->tightenPrimalBounds();
00279           if (numberInfeasibilities)
00280             std::cout<<"** Analysis indicates model infeasible"
00281                      <<std::endl;
00282           //if (doIdiot<0)
00283           //model2->crash(1000,1);
00284           model2->dual();
00285         } else {
00286           if (doIdiot>0) {
00287             Idiot info(*model2);
00288             info.crash(doIdiot,model2->messageHandler(),model2->messagesPointer());
00289           }
00290           model2->primal(1);
00291         }
00292         pinfo.postsolve(true);
00293 
00294         delete model2;
00295         solution.checkSolution();
00296         if (!solution.isProvenOptimal()) {
00297           printf("Resolving from postsolved model\n");
00298           
00299           int savePerturbation = solution.perturbation();
00300           solution.setPerturbation(100);
00301           solution.primal(1);
00302           solution.setPerturbation(savePerturbation);
00303           if (solution.numberIterations())
00304             printf("****** iterated %d\n",solution.numberIterations());
00305         }
00306         /*solution.checkSolution();
00307         printf("%g dual %g(%d) Primal %g(%d)\n",
00308                solution.objectiveValue(),
00309                solution.sumDualInfeasibilities(),
00310                solution.numberDualInfeasibilities(),
00311                solution.sumPrimalInfeasibilities(),
00312                solution.numberPrimalInfeasibilities());*/
00313         if (0) {
00314           ClpPresolve pinfoA;
00315           model2 = pinfoA.presolvedModel(solution,1.0e-8);
00316 
00317           printf("Resolving from presolved optimal solution\n");
00318           model2->primal(1);
00319 
00320           delete model2;
00321         }
00322       } else {
00323         if (doDual) {
00324           if (doIdiot<0)
00325             solution.crash(1000,1);
00326           solution.dual();
00327         } else {
00328           if (doIdiot>0) {
00329             Idiot info(solution);
00330             info.crash(doIdiot,solution.messageHandler(),solution.messagesPointer());
00331           }
00332           solution.primal(1);
00333         }
00334       }
00335       double time2 = CoinCpuTime()-time1;
00336       timeTaken += time2;
00337       printf("Took %g seconds\n",time2);
00338       // Test objective solution value
00339       {
00340         double soln = solution.objectiveValue();
00341         CoinRelFltEq eq(objValueTol[m]);
00342         std::cerr <<soln <<",  " <<objValue[m] <<" diff "<<
00343           soln-objValue[m]<<std::endl;
00344         if(!eq(soln,objValue[m]))
00345           printf("** difference fails\n");
00346       }
00347     }
00348     printf("Total time %g seconds\n",timeTaken);
00349   }
00350   else {
00351     testingMessage( "***Skipped Testing on netlib    ***\n" );
00352     testingMessage( "***use -netlib to test class***\n" );
00353   }
00354   
00355   testingMessage( "All tests completed successfully\n" );
00356   return 0;
00357 }
00358 
00359  
00360 // Display message on stdout and stderr
00361 void testingMessage( const char * const msg )
00362 {
00363   std::cerr <<msg;
00364   //cout <<endl <<"*****************************************"
00365   //     <<endl <<msg <<endl;
00366 }
00367 
00368 //--------------------------------------------------------------------------
00369 // test factorization methods and simplex method and simple barrier
00370 void
00371 ClpSimplexUnitTest(const std::string & mpsDir,
00372                    const std::string & netlibDir)
00373 {
00374   
00375   CoinRelFltEq eq(0.000001);
00376 
00377   {
00378     ClpSimplex solution;
00379   
00380     // matrix data
00381     //deliberate hiccup of 2 between 0 and 1
00382     CoinBigIndex start[5]={0,4,7,8,9};
00383     int length[5]={2,3,1,1,1};
00384     int rows[11]={0,2,-1,-1,0,1,2,0,1,2};
00385     double elements[11]={7.0,2.0,1.0e10,1.0e10,-2.0,1.0,-2.0,1,1,1};
00386     CoinPackedMatrix matrix(true,3,5,8,elements,rows,start,length);
00387     
00388     // rim data
00389     double objective[7]={-4.0,1.0,0.0,0.0,0.0,0.0,0.0};
00390     double rowLower[5]={14.0,3.0,3.0,1.0e10,1.0e10};
00391     double rowUpper[5]={14.0,3.0,3.0,-1.0e10,-1.0e10};
00392     double colLower[7]={0.0,0.0,0.0,0.0,0.0,0.0,0.0};
00393     double colUpper[7]={100.0,100.0,100.0,100.0,100.0,100.0,100.0};
00394     
00395     // basis 1
00396     int rowBasis1[3]={-1,-1,-1};
00397     int colBasis1[5]={1,1,-1,-1,1};
00398     solution.loadProblem(matrix,colLower,colUpper,objective,
00399                          rowLower,rowUpper);
00400     int i;
00401     solution.createStatus();
00402     for (i=0;i<3;i++) {
00403       if (rowBasis1[i]<0) {
00404         solution.setRowStatus(i,ClpSimplex::atLowerBound);
00405       } else {
00406         solution.setRowStatus(i,ClpSimplex::basic);
00407       }
00408     }
00409     for (i=0;i<5;i++) {
00410       if (colBasis1[i]<0) {
00411         solution.setColumnStatus(i,ClpSimplex::atLowerBound);
00412       } else {
00413         solution.setColumnStatus(i,ClpSimplex::basic);
00414       }
00415     }
00416     solution.setLogLevel(3+4+8+16+32);
00417     solution.primal();
00418     for (i=0;i<3;i++) {
00419       if (rowBasis1[i]<0) {
00420         solution.setRowStatus(i,ClpSimplex::atLowerBound);
00421       } else {
00422         solution.setRowStatus(i,ClpSimplex::basic);
00423       }
00424     }
00425     for (i=0;i<5;i++) {
00426       if (colBasis1[i]<0) {
00427         solution.setColumnStatus(i,ClpSimplex::atLowerBound);
00428       } else {
00429         solution.setColumnStatus(i,ClpSimplex::basic);
00430       }
00431     }
00432     // intricate stuff does not work with scaling
00433     solution.scaling(0);
00434     assert(!solution.factorize ( ));
00435     const double * colsol = solution.primalColumnSolution();
00436     const double * rowsol = solution.primalRowSolution();
00437     solution.getSolution(rowsol,colsol);
00438     double colsol1[5]={20.0/7.0,3.0,0.0,0.0,23.0/7.0};
00439     for (i=0;i<5;i++) {
00440       assert(eq(colsol[i],colsol1[i]));
00441     }
00442     // now feed in again without actually doing factorization
00443     ClpFactorization factorization2 = *solution.factorization();
00444     ClpSimplex solution2 = solution;
00445     solution2.setFactorization(factorization2);
00446     solution2.createStatus();
00447     for (i=0;i<3;i++) {
00448       if (rowBasis1[i]<0) {
00449         solution2.setRowStatus(i,ClpSimplex::atLowerBound);
00450       } else {
00451         solution2.setRowStatus(i,ClpSimplex::basic);
00452       }
00453     }
00454     for (i=0;i<5;i++) {
00455       if (colBasis1[i]<0) {
00456         solution2.setColumnStatus(i,ClpSimplex::atLowerBound);
00457       } else {
00458         solution2.setColumnStatus(i,ClpSimplex::basic);
00459       }
00460     }
00461     // intricate stuff does not work with scaling
00462     solution2.scaling(0);
00463     solution2.getSolution(rowsol,colsol);
00464     colsol = solution2.primalColumnSolution();
00465     rowsol = solution2.primalRowSolution();
00466     for (i=0;i<5;i++) {
00467       assert(eq(colsol[i],colsol1[i]));
00468     }
00469     solution2.setDualBound(0.1);
00470     solution2.dual();
00471     objective[2]=-1.0;
00472     objective[3]=-0.5;
00473     objective[4]=10.0;
00474     solution.dual();
00475     for (i=0;i<3;i++) {
00476       rowLower[i]=-1.0e50;
00477       colUpper[i+2]=0.0;
00478     }
00479     solution.setLogLevel(3);
00480     solution.dual();
00481     double rowObjective[]={1.0,0.5,-10.0};
00482     solution.loadProblem(matrix,colLower,colUpper,objective,
00483                          rowLower,rowUpper,rowObjective);
00484     solution.dual();
00485   }
00486   {    
00487     CoinMpsIO m;
00488     std::string fn = mpsDir+"exmip1";
00489     m.readMps(fn.c_str(),"mps");
00490     ClpSimplex solution;
00491     solution.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
00492                          m.getObjCoefficients(),
00493                          m.getRowLower(),m.getRowUpper());
00494     solution.dual();
00495   }
00496   // Test Message handler
00497   {    
00498     CoinMpsIO m;
00499     std::string fn = mpsDir+"exmip1";
00500     //fn = "Test/subGams4";
00501     m.readMps(fn.c_str(),"mps");
00502     ClpSimplex model;
00503     model.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
00504                          m.getObjCoefficients(),
00505                          m.getRowLower(),m.getRowUpper());
00506     // Message handler
00507     MyMessageHandler messageHandler(&model);
00508     std::cout<<"Testing derived message handler"<<std::endl;
00509     model.passInMessageHandler(&messageHandler);
00510     model.messagesPointer()->setDetailMessage(1,102);
00511     model.setFactorizationFrequency(10);
00512     model.primal();
00513 
00514     // Write saved solutions
00515     int nc = model.getNumCols();
00516     int s; 
00517     std::deque<StdVectorDouble> fep = messageHandler.getFeasibleExtremePoints();
00518     int numSavedSolutions = fep.size();
00519     for ( s=0; s<numSavedSolutions; ++s ) {
00520       const StdVectorDouble & solnVec = fep[s];
00521       for ( int c=0; c<nc; ++c ) {
00522         if (fabs(solnVec[c])>1.0e-8)
00523           std::cout <<"Saved Solution: " <<s <<" ColNum: " <<c <<" Value: " <<solnVec[c] <<std::endl;
00524       }
00525     }
00526     // Solve again without scaling
00527     // and maximize then minimize
00528     messageHandler.clearFeasibleExtremePoints();
00529     model.scaling(0);
00530     model.setOptimizationDirection(-1);
00531     model.primal();
00532     model.setOptimizationDirection(1);
00533     model.primal();
00534     fep = messageHandler.getFeasibleExtremePoints();
00535     numSavedSolutions = fep.size();
00536     for ( s=0; s<numSavedSolutions; ++s ) {
00537       const StdVectorDouble & solnVec = fep[s];
00538       for ( int c=0; c<nc; ++c ) {
00539         if (fabs(solnVec[c])>1.0e-8)
00540           std::cout <<"Saved Solution: " <<s <<" ColNum: " <<c <<" Value: " <<solnVec[c] <<std::endl;
00541       }
00542     }
00543   }
00544   // test steepest edge
00545   {    
00546     CoinMpsIO m;
00547     std::string fn = netlibDir+"finnis";
00548     m.readMps(fn.c_str(),"mps");
00549     ClpModel model;
00550     model.loadProblem(*m.getMatrixByCol(),m.getColLower(),
00551                     m.getColUpper(),
00552                     m.getObjCoefficients(),
00553                     m.getRowLower(),m.getRowUpper());
00554     ClpSimplex solution(model);
00555 
00556     solution.scaling(1); 
00557     solution.setDualBound(1.0e8);
00558     //solution.factorization()->maximumPivots(1);
00559     //solution.setLogLevel(3);
00560     solution.setDualTolerance(1.0e-7);
00561     // set objective sense,
00562     ClpDualRowSteepest steep;
00563     solution.setDualRowPivotAlgorithm(steep);
00564     solution.setDblParam(ClpObjOffset,m.objectiveOffset());
00565     solution.dual();
00566   }
00567   // test normal solution
00568   {    
00569     CoinMpsIO m;
00570     std::string fn = netlibDir+"afiro";
00571     m.readMps(fn.c_str(),"mps");
00572     ClpSimplex solution;
00573     ClpModel model;
00574     // do twice - without and with scaling
00575     int iPass;
00576     for (iPass=0;iPass<2;iPass++) {
00577       // explicit row objective for testing
00578       int nr = m.getNumRows();
00579       double * rowObj = new double[nr];
00580       CoinFillN(rowObj,nr,0.0);
00581       model.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
00582                       m.getObjCoefficients(),
00583                       m.getRowLower(),m.getRowUpper(),rowObj);
00584       delete [] rowObj;
00585       solution = ClpSimplex(model);
00586       if (iPass) {
00587         solution.scaling();
00588       }
00589       solution.dual();
00590       solution.dual();
00591       // test optimal
00592       assert (solution.status()==0);
00593       int numberColumns = solution.numberColumns();
00594       int numberRows = solution.numberRows();
00595       CoinPackedVector colsol(numberColumns,solution.primalColumnSolution());
00596       double * objective = solution.objective();
00597       double objValue = colsol.dotProduct(objective);
00598       CoinRelFltEq eq(1.0e-8);
00599       assert(eq(objValue,-4.6475314286e+02));
00600       double * lower = solution.columnLower();
00601       double * upper = solution.columnUpper();
00602       double * sol = solution.primalColumnSolution();
00603       double * result = new double[numberColumns];
00604       CoinFillN ( result, numberColumns,0.0);
00605       solution.matrix()->transposeTimes(solution.dualRowSolution(), result);
00606       int iRow , iColumn;
00607       // see if feasible and dual feasible
00608       for (iColumn=0;iColumn<numberColumns;iColumn++) {
00609         double value = sol[iColumn];
00610         assert(value<upper[iColumn]+1.0e-8);
00611         assert(value>lower[iColumn]-1.0e-8);
00612         value = objective[iColumn]-result[iColumn];
00613         assert (value>-1.0e-5);
00614         if (sol[iColumn]>1.0e-5)
00615           assert (value<1.0e-5);
00616       }
00617       delete [] result;
00618       result = new double[numberRows];
00619       CoinFillN ( result, numberRows,0.0);
00620       solution.matrix()->times(colsol, result);
00621       lower = solution.rowLower();
00622       upper = solution.rowUpper();
00623       sol = solution.primalRowSolution();
00624       for (iRow=0;iRow<numberRows;iRow++) {
00625         double value = result[iRow];
00626         assert(eq(value,sol[iRow]));
00627         assert(value<upper[iRow]+1.0e-8);
00628         assert(value>lower[iRow]-1.0e-8);
00629       }
00630       delete [] result;
00631       // test row objective
00632       double * rowObjective = solution.rowObjective();
00633       CoinDisjointCopyN(solution.dualRowSolution(),numberRows,rowObjective);
00634       CoinDisjointCopyN(solution.dualColumnSolution(),numberColumns,objective);
00635       // this sets up all slack basis
00636       solution.createStatus();
00637       solution.dual();
00638       CoinFillN(rowObjective,numberRows,0.0);
00639       CoinDisjointCopyN(m.getObjCoefficients(),numberColumns,objective);
00640       solution.dual();
00641     }
00642   }
00643   // test unbounded
00644   {    
00645     CoinMpsIO m;
00646     std::string fn = netlibDir+"brandy";
00647     m.readMps(fn.c_str(),"mps");
00648     ClpSimplex solution;
00649     // do twice - without and with scaling
00650     int iPass;
00651     for (iPass=0;iPass<2;iPass++) {
00652       solution.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
00653                       m.getObjCoefficients(),
00654                       m.getRowLower(),m.getRowUpper());
00655       if (iPass)
00656         solution.scaling();
00657       solution.setOptimizationDirection(-1);
00658       // test unbounded and ray
00659 #ifdef DUAL
00660       solution.setDualBound(100.0);
00661       solution.dual();
00662 #else
00663       solution.primal();
00664 #endif
00665       assert (solution.status()==2);
00666       int numberColumns = solution.numberColumns();
00667       int numberRows = solution.numberRows();
00668       double * lower = solution.columnLower();
00669       double * upper = solution.columnUpper();
00670       double * sol = solution.primalColumnSolution();
00671       double * ray = solution.unboundedRay();
00672       double * objective = solution.objective();
00673       double objChange=0.0;
00674       int iRow , iColumn;
00675       // make sure feasible and columns form ray
00676       for (iColumn=0;iColumn<numberColumns;iColumn++) {
00677         double value = sol[iColumn];
00678         assert(value<upper[iColumn]+1.0e-8);
00679         assert(value>lower[iColumn]-1.0e-8);
00680         value = ray[iColumn];
00681         if (value>0.0)
00682           assert(upper[iColumn]>1.0e30);
00683         else if (value<0.0)
00684           assert(lower[iColumn]<-1.0e30);
00685         objChange += value*objective[iColumn];
00686       }
00687       // make sure increasing objective
00688       assert(objChange>0.0);
00689       double * result = new double[numberRows];
00690       CoinFillN ( result, numberRows,0.0);
00691       solution.matrix()->times(sol, result);
00692       lower = solution.rowLower();
00693       upper = solution.rowUpper();
00694       sol = solution.primalRowSolution();
00695       for (iRow=0;iRow<numberRows;iRow++) {
00696         double value = result[iRow];
00697         assert(eq(value,sol[iRow]));
00698         assert(value<upper[iRow]+2.0e-8);
00699         assert(value>lower[iRow]-2.0e-8);
00700       }
00701       CoinFillN ( result, numberRows,0.0);
00702       solution.matrix()->times(ray, result);
00703       // there may be small differences (especially if scaled)
00704       for (iRow=0;iRow<numberRows;iRow++) {
00705         double value = result[iRow];
00706         if (value>1.0e-8)
00707           assert(upper[iRow]>1.0e30);
00708         else if (value<-1.0e-8)
00709           assert(lower[iRow]<-1.0e30);
00710       }
00711       delete [] result;
00712       delete [] ray;
00713     }
00714   }
00715   // test infeasible
00716   {    
00717     CoinMpsIO m;
00718     std::string fn = netlibDir+"brandy";
00719     m.readMps(fn.c_str(),"mps");
00720     ClpSimplex solution;
00721     // do twice - without and with scaling
00722     int iPass;
00723     for (iPass=0;iPass<2;iPass++) {
00724       solution.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
00725                       m.getObjCoefficients(),
00726                       m.getRowLower(),m.getRowUpper());
00727       if (iPass)
00728         solution.scaling();
00729       // test infeasible and ray
00730       solution.columnUpper()[0]=0.0;
00731 #ifdef DUAL
00732       solution.setDualBound(100.0);
00733       solution.dual();
00734 #else
00735       solution.primal();
00736 #endif
00737       assert (solution.status()==1);
00738       int numberColumns = solution.numberColumns();
00739       int numberRows = solution.numberRows();
00740       double * lower = solution.rowLower();
00741       double * upper = solution.rowUpper();
00742       double * ray = solution.infeasibilityRay();
00743       assert(ray);
00744       // construct proof of infeasibility
00745       int iRow , iColumn;
00746       double lo=0.0,up=0.0;
00747       int nl=0,nu=0;
00748       for (iRow=0;iRow<numberRows;iRow++) {
00749         if (lower[iRow]>-1.0e20) {
00750           lo += ray[iRow]*lower[iRow];
00751         } else {
00752           if (ray[iRow]>1.0e-8) 
00753             nl++;
00754         }
00755         if (upper[iRow]<1.0e20) {
00756           up += ray[iRow]*upper[iRow];
00757         } else {
00758           if (ray[iRow]>1.0e-8) 
00759             nu++;
00760         }
00761       }
00762       if (nl)
00763         lo=-1.0e100;
00764       if (nu)
00765         up=1.0e100;
00766       double * result = new double[numberColumns];
00767       double lo2=0.0,up2=0.0;
00768       CoinFillN ( result, numberColumns,0.0);
00769       solution.matrix()->transposeTimes(ray, result);
00770       lower = solution.columnLower();
00771       upper = solution.columnUpper();
00772       nl=nu=0;
00773       for (iColumn=0;iColumn<numberColumns;iColumn++) {
00774         if (result[iColumn]>1.0e-8) {
00775           if (lower[iColumn]>-1.0e20)
00776             lo2 += result[iColumn]*lower[iColumn];
00777           else
00778             nl++;
00779           if (upper[iColumn]<1.0e20)
00780             up2 += result[iColumn]*upper[iColumn];
00781           else
00782             nu++;
00783         } else if (result[iColumn]<-1.0e-8) {
00784           if (lower[iColumn]>-1.0e20)
00785             up2 += result[iColumn]*lower[iColumn];
00786           else
00787             nu++;
00788           if (upper[iColumn]<1.0e20)
00789             lo2 += result[iColumn]*upper[iColumn];
00790           else
00791             nl++;
00792         }
00793       }
00794       if (nl)
00795         lo2=-1.0e100;
00796       if (nu)
00797         up2=1.0e100;
00798       // make sure inconsistency
00799       assert(lo2>up||up2<lo);
00800       delete [] result;
00801       delete [] ray;
00802     }
00803   }
00804   // test delete and add
00805   {    
00806     CoinMpsIO m;
00807     std::string fn = netlibDir+"brandy";
00808     m.readMps(fn.c_str(),"mps");
00809     ClpSimplex solution;
00810     solution.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
00811                          m.getObjCoefficients(),
00812                          m.getRowLower(),m.getRowUpper());
00813     solution.dual();
00814     CoinRelFltEq eq(1.0e-8);
00815     assert(eq(solution.objectiveValue(),1.5185098965e+03));
00816 
00817     int numberColumns = solution.numberColumns();
00818     int numberRows = solution.numberRows();
00819     double * saveObj = new double [numberColumns];
00820     double * saveLower = new double[numberRows+numberColumns];
00821     double * saveUpper = new double[numberRows+numberColumns];
00822     int * which = new int [numberRows+numberColumns];
00823 
00824     int numberElements = m.getMatrixByCol()->getNumElements();
00825     int * starts = new int[numberRows+numberColumns];
00826     int * index = new int[numberElements];
00827     double * element = new double[numberElements];
00828 
00829     const int * startM;
00830     const int * lengthM;
00831     const int * indexM;
00832     const double * elementM;
00833 
00834     int n,nel;
00835 
00836     // delete non basic columns
00837     n=0;
00838     nel=0;
00839     int iRow , iColumn;
00840     double * lower = solution.columnLower();
00841     double * upper = solution.columnUpper();
00842     double * objective = solution.objective();
00843     startM = m.getMatrixByCol()->getVectorStarts();
00844     lengthM = m.getMatrixByCol()->getVectorLengths();
00845     indexM = m.getMatrixByCol()->getIndices();
00846     elementM = m.getMatrixByCol()->getElements();
00847     starts[0]=0;
00848     for (iColumn=0;iColumn<numberColumns;iColumn++) {
00849       if (solution.getColumnStatus(iColumn)!=ClpSimplex::basic) {
00850         saveObj[n]=objective[iColumn];
00851         saveLower[n]=lower[iColumn];
00852         saveUpper[n]=upper[iColumn];
00853         int j;
00854         for (j=startM[iColumn];j<startM[iColumn]+lengthM[iColumn];j++) {
00855           index[nel]=indexM[j];
00856           element[nel++]=elementM[j];
00857         }
00858         which[n++]=iColumn;
00859         starts[n]=nel;
00860       }
00861     }
00862     solution.deleteColumns(n,which);
00863     solution.dual();
00864     // Put back
00865     solution.addColumns(n,saveLower,saveUpper,saveObj,
00866                         starts,index,element);
00867     solution.dual();
00868     assert(eq(solution.objectiveValue(),1.5185098965e+03));
00869 
00870     // reload with original
00871     solution.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
00872                          m.getObjCoefficients(),
00873                          m.getRowLower(),m.getRowUpper());
00874     // delete half rows
00875     n=0;
00876     nel=0;
00877     lower = solution.rowLower();
00878     upper = solution.rowUpper();
00879     startM = m.getMatrixByRow()->getVectorStarts();
00880     lengthM = m.getMatrixByRow()->getVectorLengths();
00881     indexM = m.getMatrixByRow()->getIndices();
00882     elementM = m.getMatrixByRow()->getElements();
00883     starts[0]=0;
00884     for (iRow=0;iRow<numberRows;iRow++) {
00885       if ((iRow&1)==0) {
00886         saveLower[n]=lower[iRow];
00887         saveUpper[n]=upper[iRow];
00888         int j;
00889         for (j=startM[iRow];j<startM[iRow]+lengthM[iRow];j++) {
00890           index[nel]=indexM[j];
00891           element[nel++]=elementM[j];
00892         }
00893         which[n++]=iRow;
00894         starts[n]=nel;
00895       }
00896     }
00897     solution.deleteRows(n,which);
00898     solution.dual();
00899     // Put back
00900     solution.addRows(n,saveLower,saveUpper,
00901                         starts,index,element);
00902     solution.dual();
00903     assert(eq(solution.objectiveValue(),1.5185098965e+03));
00904     // Zero out status array to give some interest
00905     memset(solution.statusArray()+numberColumns,0,numberRows);
00906     solution.primal(1);
00907     assert(eq(solution.objectiveValue(),1.5185098965e+03));
00908 
00909     delete [] saveObj;
00910     delete [] saveLower;
00911     delete [] saveUpper;
00912     delete [] which;
00913     delete [] starts;
00914     delete [] index;
00915     delete [] element;
00916   }
00917 #if 0
00918   // Test barrier
00919   {
00920     CoinMpsIO m;
00921     std::string fn = mpsDir+"exmip1";
00922     m.readMps(fn.c_str(),"mps");
00923     ClpInterior solution;
00924     solution.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
00925                          m.getObjCoefficients(),
00926                          m.getRowLower(),m.getRowUpper());
00927     solution.primalDual();
00928   }
00929 #endif
00930   // test network 
00931   //#define QUADRATIC
00932 #ifndef QUADRATIC
00933   if (1) {    
00934     std::string fn = mpsDir+"input.130";
00935     int numberColumns;
00936     int numberRows;
00937     
00938     FILE * fp = fopen(fn.c_str(),"r");
00939     if (!fp) {
00940       // Try in Samples
00941       fn = "Samples/input.130";
00942       fp = fopen(fn.c_str(),"r");
00943     }
00944     if (!fp) {
00945       fprintf(stderr,"Unable to open file input.130 in mpsDir or Samples directory\n");
00946     } else {
00947       int problem;
00948       char temp[100];
00949       // read and skip 
00950       fscanf(fp,"%s",temp);
00951       assert (!strcmp(temp,"BEGIN"));
00952       fscanf(fp,"%*s %*s %d %d %*s %*s %d %*s",&problem, &numberRows, 
00953              &numberColumns);
00954       // scan down to SUPPLY
00955       while (fgets(temp,100,fp)) {
00956         if (!strncmp(temp,"SUPPLY",6))
00957           break;
00958       }
00959       if (strncmp(temp,"SUPPLY",6)) {
00960         fprintf(stderr,"Unable to find SUPPLY\n");
00961         exit(2);
00962       }
00963       // get space for rhs
00964       double * lower = new double[numberRows];
00965       double * upper = new double[numberRows];
00966       int i;
00967       for (i=0;i<numberRows;i++) {
00968         lower[i]=0.0;
00969         upper[i]=0.0;
00970       }
00971       // ***** Remember to convert to C notation
00972       while (fgets(temp,100,fp)) {
00973         int row;
00974         int value;
00975         if (!strncmp(temp,"ARCS",4))
00976           break;
00977         sscanf(temp,"%d %d",&row,&value);
00978         upper[row-1]=-value;
00979         lower[row-1]=-value;
00980       }
00981       if (strncmp(temp,"ARCS",4)) {
00982         fprintf(stderr,"Unable to find ARCS\n");
00983         exit(2);
00984       }
00985       // number of columns may be underestimate
00986       int * head = new int[2*numberColumns];
00987       int * tail = new int[2*numberColumns];
00988       int * ub = new int[2*numberColumns];
00989       int * cost = new int[2*numberColumns];
00990       // ***** Remember to convert to C notation
00991       numberColumns=0;
00992       while (fgets(temp,100,fp)) {
00993         int iHead;
00994         int iTail;
00995         int iUb;
00996         int iCost;
00997         if (!strncmp(temp,"DEMAND",6))
00998           break;
00999         sscanf(temp,"%d %d %d %d",&iHead,&iTail,&iCost,&iUb);
01000         iHead--;
01001         iTail--;
01002         head[numberColumns]=iHead;
01003         tail[numberColumns]=iTail;
01004         ub[numberColumns]=iUb;
01005         cost[numberColumns]=iCost;
01006         numberColumns++;
01007       }
01008       if (strncmp(temp,"DEMAND",6)) {
01009         fprintf(stderr,"Unable to find DEMAND\n");
01010         exit(2);
01011       }
01012       // ***** Remember to convert to C notation
01013       while (fgets(temp,100,fp)) {
01014         int row;
01015         int value;
01016         if (!strncmp(temp,"END",3))
01017           break;
01018         sscanf(temp,"%d %d",&row,&value);
01019         upper[row-1]=value;
01020         lower[row-1]=value;
01021       }
01022       if (strncmp(temp,"END",3)) {
01023         fprintf(stderr,"Unable to find END\n");
01024         exit(2);
01025       }
01026       printf("Problem %d has %d rows and %d columns\n",problem,
01027              numberRows,numberColumns);
01028       fclose(fp);
01029       ClpSimplex  model;
01030       // now build model
01031       
01032       double * objective =new double[numberColumns];
01033       double * lowerColumn = new double[numberColumns];
01034       double * upperColumn = new double[numberColumns];
01035       
01036       double * element = new double [2*numberColumns];
01037       int * start = new int[numberColumns+1];
01038       int * row = new int[2*numberColumns];
01039       start[numberColumns]=2*numberColumns;
01040       for (i=0;i<numberColumns;i++) {
01041         start[i]=2*i;
01042         element[2*i]=-1.0;
01043         element[2*i+1]=1.0;
01044         row[2*i]=head[i];
01045         row[2*i+1]=tail[i];
01046         lowerColumn[i]=0.0;
01047         upperColumn[i]=ub[i];
01048         objective[i]=cost[i];
01049       }
01050       // Create Packed Matrix
01051       CoinPackedMatrix matrix;
01052       int * lengths = NULL;
01053       matrix.assignMatrix(true,numberRows,numberColumns,
01054                           2*numberColumns,element,row,start,lengths);
01055       // load model
01056       
01057       model.loadProblem(matrix,
01058                         lowerColumn,upperColumn,objective,
01059                         lower,upper);
01060       
01061       model.factorization()->maximumPivots(200+model.numberRows()/100);
01062       model.createStatus();
01063       double time1 = CoinCpuTime();
01064       model.dual();
01065       std::cout<<"Network problem, ClpPackedMatrix took "<<CoinCpuTime()-time1<<" seconds"<<std::endl;
01066       ClpPlusMinusOneMatrix plusMinus(matrix);
01067       assert (plusMinus.getIndices()); // would be zero if not +- one
01068       model.loadProblem(plusMinus,
01069                         lowerColumn,upperColumn,objective,
01070                         lower,upper);
01071       
01072       model.factorization()->maximumPivots(200+model.numberRows()/100);
01073       model.createStatus();
01074       time1 = CoinCpuTime();
01075       model.dual();
01076       std::cout<<"Network problem, ClpPlusMinusOneMatrix took "<<CoinCpuTime()-time1<<" seconds"<<std::endl;
01077       ClpNetworkMatrix network(numberColumns,head,tail);
01078       model.loadProblem(network,
01079                         lowerColumn,upperColumn,objective,
01080                         lower,upper);
01081       
01082       model.factorization()->maximumPivots(200+model.numberRows()/100);
01083       model.createStatus();
01084       time1 = CoinCpuTime();
01085       model.dual();
01086       std::cout<<"Network problem, ClpNetworkMatrix took "<<CoinCpuTime()-time1<<" seconds"<<std::endl;
01087       delete [] lower;
01088       delete [] upper;
01089       delete [] head;
01090       delete [] tail;
01091       delete [] ub;
01092       delete [] cost;
01093       delete [] objective;
01094       delete [] lowerColumn;
01095       delete [] upperColumn;
01096     }
01097   }
01098 #endif
01099 #ifdef QUADRATIC
01100   // Test quadratic to solve linear
01101   if (1) {    
01102     CoinMpsIO m;
01103     std::string fn = mpsDir+"exmip1";
01104     m.readMps(fn.c_str(),"mps");
01105     ClpSimplex solution;
01106     solution.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
01107                          m.getObjCoefficients(),
01108                          m.getRowLower(),m.getRowUpper());
01109     //solution.dual();
01110     // get quadratic part
01111     int numberColumns=solution.numberColumns();
01112     int * start=new int [numberColumns+1];
01113     int * column = new int[numberColumns];
01114     double * element = new double[numberColumns];
01115     int i;
01116     start[0]=0;
01117     int n=0;
01118     int kk=numberColumns-1;
01119     int kk2=numberColumns-1;
01120     for (i=0;i<numberColumns;i++) {
01121       if (i>=kk) {
01122         column[n]=i;
01123         if (i>=kk2)
01124           element[n]=1.0e-1;
01125         else
01126           element[n]=0.0;
01127         n++;
01128       }
01129       start[i+1]=n;
01130     }
01131     // Load up objective
01132     solution.loadQuadraticObjective(numberColumns,start,column,element);
01133     delete [] start;
01134     delete [] column;
01135     delete [] element;
01136     //solution.quadraticSLP(50,1.0e-4);
01137     double objValue = solution.getObjValue();
01138     CoinRelFltEq eq(1.0e-4);
01139     //assert(eq(objValue,820.0));
01140     solution.setLogLevel(63);
01141     solution.quadraticPrimal(1);
01142     objValue = solution.getObjValue();
01143     //assert(eq(objValue,3.2368421));
01144     //exit(77);
01145   }
01146   // Test quadratic
01147   if (1) {    
01148     CoinMpsIO m;
01149     std::string fn = mpsDir+"share2qp";
01150     //fn = "share2qpb";
01151     m.readMps(fn.c_str(),"mps");
01152     ClpSimplex model;
01153     model.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
01154                          m.getObjCoefficients(),
01155                          m.getRowLower(),m.getRowUpper());
01156     model.dual();
01157     // get quadratic part
01158     int * start=NULL;
01159     int * column = NULL;
01160     double * element = NULL;
01161     m.readQuadraticMps(NULL,start,column,element,2);
01162     int column2[200];
01163     double element2[200];
01164     int start2[80];
01165     int j;
01166     start2[0]=0;
01167     int nel=0;
01168     bool good=false;
01169     for (j=0;j<79;j++) {
01170       if (start[j]==start[j+1]) {
01171         column2[nel]=j;
01172         element2[nel]=0.0;
01173         nel++;
01174       } else {
01175         int i;
01176         for (i=start[j];i<start[j+1];i++) {
01177           column2[nel]=column[i];
01178           element2[nel++]=element[i];
01179         }
01180       }
01181       start2[j+1]=nel;
01182     }
01183     // Load up objective
01184     if (good)
01185       model.loadQuadraticObjective(model.numberColumns(),start2,column2,element2);
01186     else
01187       model.loadQuadraticObjective(model.numberColumns(),start,column,element);
01188     delete [] start;
01189     delete [] column;
01190     delete [] element;
01191     int numberColumns=model.numberColumns();
01192 #if 0
01193     model.quadraticSLP(50,1.0e-4);
01194 #else
01195     // Get feasible
01196     ClpObjective * saveObjective = model.objectiveAsObject()->clone();
01197     ClpLinearObjective zeroObjective(NULL,numberColumns);
01198     model.setObjective(&zeroObjective);
01199     model.dual();
01200     model.setObjective(saveObjective);
01201     delete saveObjective;
01202 #endif
01203     model.setLogLevel(63);
01204     //exit(77);
01205     model.setFactorizationFrequency(10);
01206     model.quadraticPrimal(1);
01207     double objValue = model.getObjValue();
01208     const double * solution = model.primalColumnSolution();
01209     int i;
01210     for (i=0;i<numberColumns;i++)
01211       if (solution[i])
01212         printf("%d %g\n",i,solution[i]);
01213     CoinRelFltEq eq(1.0e-4);
01214     assert(eq(objValue,-400.92));
01215   }
01216   if (0) {    
01217     CoinMpsIO m;
01218     std::string fn = "./beale";
01219     //fn = "./jensen";
01220     m.readMps(fn.c_str(),"mps");
01221     ClpSimplex solution;
01222     solution.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
01223                          m.getObjCoefficients(),
01224                          m.getRowLower(),m.getRowUpper());
01225     solution.setDblParam(ClpObjOffset,m.objectiveOffset());
01226     solution.dual();
01227     // get quadratic part
01228     int * start=NULL;
01229     int * column = NULL;
01230     double * element = NULL;
01231     m.readQuadraticMps(NULL,start,column,element,2);
01232     // Load up objective
01233     solution.loadQuadraticObjective(solution.numberColumns(),start,column,element);
01234     delete [] start;
01235     delete [] column;
01236     delete [] element;
01237     solution.quadraticPrimal(1);
01238     solution.quadraticSLP(50,1.0e-4);
01239     double objValue = solution.getObjValue();
01240     CoinRelFltEq eq(1.0e-4);
01241     assert(eq(objValue,0.5));
01242     solution.quadraticPrimal();
01243     objValue = solution.getObjValue();
01244     assert(eq(objValue,0.5));
01245   }
01246 #endif  
01247 }