// Copyright (C) 2002, International Business Machines
// Corporation and others.  All Rights Reserved.

#include "ClpConfig.h"
#include "CoinPragma.hpp"
#include <cassert>
#include <cstdio>
#include <cmath>
#include <cfloat>
#include <string>
#include <iostream>

#include "CoinMpsIO.hpp"
#include "CoinPackedMatrix.hpp"
#include "CoinPackedVector.hpp"
#include "CoinHelperFunctions.hpp"
#include "CoinTime.hpp"

#include "ClpFactorization.hpp"
#include "ClpSimplex.hpp"
#include "ClpSimplexOther.hpp"
#include "ClpSimplexNonlinear.hpp"
#include "ClpInterior.hpp"
#include "ClpLinearObjective.hpp"
#include "ClpDualRowSteepest.hpp"
#include "ClpDualRowDantzig.hpp"
#include "ClpPrimalColumnSteepest.hpp"
#include "ClpPrimalColumnDantzig.hpp"
#include "ClpParameters.hpp"
#include "ClpNetworkMatrix.hpp"
#include "ClpPlusMinusOneMatrix.hpp"
#include "MyMessageHandler.hpp"
#include "MyEventHandler.hpp"

#include "ClpPresolve.hpp"
#include "Idiot.hpp"


//#############################################################################

#ifdef NDEBUG
#undef NDEBUG
#endif

// Function Prototypes. Function definitions is in this file.
void testingMessage( const char * const msg );
#if UFL_BARRIER
static int barrierAvailable=1;
static std::string nameBarrier="barrier-UFL";
#elif WSSMP_BARRIER
static int barrierAvailable=2;
static std::string nameBarrier="barrier-WSSMP";
#else
static int barrierAvailable=0;
static std::string nameBarrier="barrier-slow";
#endif
#define NUMBER_ALGORITHMS 12
// If you just want a subset then set some to 1
static int switchOff[NUMBER_ALGORITHMS]={0,0,0,0,0,0,0,0,0,0,0,0};
// shortName - 0 no , 1 yes
ClpSolve setupForSolve(int algorithm, std::string & nameAlgorithm,
                       int shortName)
{
  ClpSolve solveOptions;
  /* algorithms are
     0 barrier
     1 dual with volumne crash
     2,3 dual with and without crash
     4,5 primal with and without
     6,7 automatic with and without
     8,9 primal with idiot 1 and 5
     10,11 primal with 70, dual with volume
  */
  switch (algorithm) {
  case 0:
    if (shortName)
      nameAlgorithm="ba";
    else
      nameAlgorithm="nameBarrier";
    solveOptions.setSolveType(ClpSolve::useBarrier);
    if (barrierAvailable==1) 
      solveOptions.setSpecialOption(4,4);
    else if (barrierAvailable==2) 
      solveOptions.setSpecialOption(4,2);
    break;
  case 1:
#ifdef COIN_HAS_VOL
    if (shortName)
      nameAlgorithm="du-vol-50";
    else
      nameAlgorithm="dual-volume-50";
    solveOptions.setSolveType(ClpSolve::useDual);
    solveOptions.setSpecialOption(0,2,50); // volume
#else  
      solveOptions.setSolveType(ClpSolve::notImplemented);
#endif
    break;
  case 2:
    if (shortName)
      nameAlgorithm="du-cr";
    else
      nameAlgorithm="dual-crash";
    solveOptions.setSolveType(ClpSolve::useDual);
    solveOptions.setSpecialOption(0,1);
    break;
  case 3:
    if (shortName)
      nameAlgorithm="du";
    else
      nameAlgorithm="dual";
    solveOptions.setSolveType(ClpSolve::useDual);
    break;
  case 4:
    if (shortName)
      nameAlgorithm="pr-cr";
    else
      nameAlgorithm="primal-crash";
    solveOptions.setSolveType(ClpSolve::usePrimal);
    solveOptions.setSpecialOption(1,1);
    break;
  case 5:
    if (shortName)
      nameAlgorithm="pr";
    else
      nameAlgorithm="primal";
    solveOptions.setSolveType(ClpSolve::usePrimal);
    break;
  case 6:
    if (shortName)
      nameAlgorithm="au-cr";
    else
      nameAlgorithm="either-crash";
    solveOptions.setSolveType(ClpSolve::automatic);
    solveOptions.setSpecialOption(1,1);
    break;
  case 7:
    if (shortName)
      nameAlgorithm="au";
    else
      nameAlgorithm="either";
    solveOptions.setSolveType(ClpSolve::automatic);
    break;
  case 8:
    if (shortName)
      nameAlgorithm="pr-id-1";
    else
      nameAlgorithm="primal-idiot-1";
    solveOptions.setSolveType(ClpSolve::usePrimalorSprint);
    solveOptions.setSpecialOption(1,2,1); // idiot
    break;
  case 9:
    if (shortName)
      nameAlgorithm="pr-id-5";
    else
      nameAlgorithm="primal-idiot-5";
    solveOptions.setSolveType(ClpSolve::usePrimalorSprint);
    solveOptions.setSpecialOption(1,2,5); // idiot
    break;
  case 10:
    if (shortName)
      nameAlgorithm="pr-id-70";
    else
      nameAlgorithm="primal-idiot-70";
    solveOptions.setSolveType(ClpSolve::usePrimalorSprint);
    solveOptions.setSpecialOption(1,2,70); // idiot
    break;
  case 11:
#ifdef COIN_HAS_VOL
    if (shortName)
      nameAlgorithm="du-vol";
    else
      nameAlgorithm="dual-volume";
    solveOptions.setSolveType(ClpSolve::useDual);
    solveOptions.setSpecialOption(0,2,3000); // volume
#else  
      solveOptions.setSolveType(ClpSolve::notImplemented);
#endif
    break;
  default:
    abort();
  }
  if (shortName) {
    // can switch off
    if (switchOff[algorithm])
      solveOptions.setSolveType(ClpSolve::notImplemented);
  }
  return solveOptions;
}
static void printSol(ClpSimplex & model) 
{
  int numberRows = model.numberRows();
  int numberColumns = model.numberColumns();
  
  double * rowPrimal = model.primalRowSolution();
  double * rowDual = model.dualRowSolution();
  double * rowLower = model.rowLower();
  double * rowUpper = model.rowUpper();
  
  int iRow;
  double objValue = model.getObjValue();
  printf("Objvalue %g Rows (%d)\n",objValue,numberRows);
  for (iRow=0;iRow<numberRows;iRow++) {
    printf("%d primal %g dual %g low %g up %g\n",
           iRow,rowPrimal[iRow],rowDual[iRow],
           rowLower[iRow],rowUpper[iRow]);
  }
  double * columnPrimal = model.primalColumnSolution();
  double * columnDual = model.dualColumnSolution();
  double * columnLower = model.columnLower();
  double * columnUpper = model.columnUpper();
  double offset;
  //const double * gradient = model.objectiveAsObject()->gradient(&model,
  //                                                       columnPrimal,offset,true,1);
  const double * gradient = model.objective(columnPrimal,offset);
  int iColumn;
  objValue = -offset-model.objectiveOffset();
  printf("offset %g (%g)\n",offset,model.objectiveOffset());
  printf("Columns (%d)\n",numberColumns);
  for (iColumn=0;iColumn<numberColumns;iColumn++) {
    printf("%d primal %g dual %g low %g up %g\n",
           iColumn,columnPrimal[iColumn],columnDual[iColumn],
           columnLower[iColumn],columnUpper[iColumn]);
    objValue += columnPrimal[iColumn]*gradient[iColumn];
    if (fabs(columnPrimal[iColumn]*gradient[iColumn])>1.0e-8)
      printf("obj -> %g gradient %g\n",objValue,gradient[iColumn]);
  }
  printf("Computed objective %g\n",objValue);
}

void usage(const std::string& key)
{
  std::cerr
    <<"Undefined parameter \"" <<key <<"\".\n"
    <<"Correct usage: \n"
    <<"  clp [-dirSample=V1] [-dirNetlib=V2] [-netlib]\n"
    <<"  where:\n"
    <<"    -dirSample: directory containing mps test files\n"
    <<"        Default value V1=\"../../Data/Sample\"\n"
    <<"    -dirNetlib: directory containing netlib files\"\n"
    <<"        Default value V2=\"../../Data/Netlib\"\n"
    <<"    -netlib\n"
    <<"        If specified, then netlib testset run as well as the nitTest.\n";
}

//----------------------------------------------------------------
int mainTest (int argc, const char *argv[],int algorithm,
	      ClpSimplex empty, bool doPresolve, int switchOffValue,bool doVector)
{
  int i;

  if (switchOffValue>0) {
    // switch off some
    int iTest;
    for (iTest=0;iTest<NUMBER_ALGORITHMS;iTest++) {
      int bottom = switchOffValue%10;
      assert (bottom==0||bottom==1);
      switchOffValue/= 10;
      switchOff[iTest]=0;
    }
  }

  // define valid parameter keywords
  std::set<std::string> definedKeyWords;
  definedKeyWords.insert("-dirSample");
  definedKeyWords.insert("-dirNetlib");
  definedKeyWords.insert("-netlib");

  // Create a map of parameter keys and associated data
  std::map<std::string,std::string> parms;
  for ( i=1; i<argc; i++ ) {
    std::string parm(argv[i]);
    std::string key,value;
    std::string::size_type  eqPos = parm.find('=');

    // Does parm contain and '='
    if ( eqPos==std::string::npos ) {
      //Parm does not contain '='
      key = parm;
    }
    else {
      key=parm.substr(0,eqPos);
      value=parm.substr(eqPos+1);
    }

    // Is specifed key valid?
    if ( definedKeyWords.find(key) == definedKeyWords.end() ) {
      // invalid key word.
      // Write help text
      usage(key);
      return 1;
    }
    parms[key]=value;
  }
  
  const char dirsep =  CoinFindDirSeparator();
  // Set directory containing mps data files.
  std::string dirSample;
  if (parms.find("-dirSample") != parms.end())
    dirSample=parms["-dirSample"];
  else 
    dirSample = dirsep == '/' ? "../../Data/Sample/" : "..\\..\\Data\\Sample\\";
 
  // Set directory containing netlib data files.
  std::string dirNetlib;
  if (parms.find("-dirNetlib") != parms.end())
    dirNetlib=parms["-dirNetlib"];
  else 
    dirNetlib = dirsep == '/' ? "../../Data/Netlib/" : "..\\..\\Data\\Netlib\\";
  if (!empty.numberRows()) {
    testingMessage( "Testing ClpSimplex\n" );
    ClpSimplexUnitTest(dirSample);
  }
  if (parms.find("-netlib") != parms.end()||empty.numberRows())
  {
    unsigned int m;
    
    // Define test problems: 
    //   mps names, 
    //   maximization or minimization, 
    //   Number of rows and columns in problem, and
    //   objective function value
    std::vector<std::string> mpsName;
    std::vector<bool> min;
    std::vector<int> nRows;
    std::vector<int> nCols;
    std::vector<double> objValue;
    std::vector<double> objValueTol;
    // 100 added means no presolve
    std::vector<int> bestStrategy;
    if(empty.numberRows()) {
      std::string alg;
      for (int iTest=0;iTest<NUMBER_ALGORITHMS;iTest++) {
        ClpSolve solveOptions=setupForSolve(iTest,alg,0);
        printf("%d %s ",iTest,alg.c_str());
        if (switchOff[iTest]) 
          printf("skipped by user\n");
        else if(solveOptions.getSolveType()==ClpSolve::notImplemented)
          printf("skipped as not available\n");
        else
          printf("will be tested\n");
      }
    }
    if (!empty.numberRows()) {
      mpsName.push_back("25fv47");
      min.push_back(true);
      nRows.push_back(822);
      nCols.push_back(1571);
      objValueTol.push_back(1.E-10);
      objValue.push_back(5.5018458883E+03);
      bestStrategy.push_back(0);
      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);bestStrategy.push_back(3);
      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);bestStrategy.push_back(3);
      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);bestStrategy.push_back(3);
      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);bestStrategy.push_back(2);
      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);bestStrategy.push_back(3);
      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);bestStrategy.push_back(3);
      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);bestStrategy.push_back(0);
      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);bestStrategy.push_back(3);
      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);bestStrategy.push_back(3);
      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);bestStrategy.push_back(3);
      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);bestStrategy.push_back(3);
      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);bestStrategy.push_back(4);
      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);bestStrategy.push_back(2);
      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);bestStrategy.push_back(0);
      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);bestStrategy.push_back(3);
      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);bestStrategy.push_back(3);
      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);bestStrategy.push_back(3);
      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);bestStrategy.push_back(3);
      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);bestStrategy.push_back(3);
      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);bestStrategy.push_back(3);
      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);bestStrategy.push_back(3);
      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);bestStrategy.push_back(3);
      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);bestStrategy.push_back(3);
      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);bestStrategy.push_back(0);
      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);bestStrategy.push_back(3);
      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);bestStrategy.push_back(3);
      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);bestStrategy.push_back(2);
      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);bestStrategy.push_back(5);
      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);bestStrategy.push_back(3); // The correct answer includes -7.113 term. This is a constant in the objective function. See line 1683 of the mps file.
      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 );bestStrategy.push_back(3);
      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);bestStrategy.push_back(3);
      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);bestStrategy.push_back(3);
      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);bestStrategy.push_back(3+100);
      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);bestStrategy.push_back(5+100);
      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);bestStrategy.push_back(3+100);
      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);bestStrategy.push_back(5+100);
      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);bestStrategy.push_back(3);
      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);bestStrategy.push_back(3);
      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);bestStrategy.push_back(3);
      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);bestStrategy.push_back(3);
      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);bestStrategy.push_back(3);
      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);bestStrategy.push_back(4+100);
      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);bestStrategy.push_back(4+100);
      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);bestStrategy.push_back(4+100);
      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);bestStrategy.push_back(2);
      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);bestStrategy.push_back(3);
      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);bestStrategy.push_back(3);
      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);bestStrategy.push_back(3);
      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);bestStrategy.push_back(3);
      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);bestStrategy.push_back(2);
      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);bestStrategy.push_back(3);
      //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);bestStrategy.push_back(3);
      //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);bestStrategy.push_back(3);
      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);bestStrategy.push_back(3);
      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);bestStrategy.push_back(3);
      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);bestStrategy.push_back(3);
      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);bestStrategy.push_back(3);
      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);bestStrategy.push_back(3);
      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);bestStrategy.push_back(3);
      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);bestStrategy.push_back(3);
      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);bestStrategy.push_back(3);
      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);bestStrategy.push_back(3);
      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);bestStrategy.push_back(3);
      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);bestStrategy.push_back(3);
      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);bestStrategy.push_back(2);
      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);bestStrategy.push_back(3+100);
      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);bestStrategy.push_back(3+100);
      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);bestStrategy.push_back(1+100);
      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);bestStrategy.push_back(3);
      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);bestStrategy.push_back(3);
      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);bestStrategy.push_back(3);
      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);bestStrategy.push_back(3);
      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);bestStrategy.push_back(3);
      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);bestStrategy.push_back(3);
      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);bestStrategy.push_back(3);
      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);bestStrategy.push_back(3);
      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);bestStrategy.push_back(3);
      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);bestStrategy.push_back(3);
      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);bestStrategy.push_back(2);
      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);bestStrategy.push_back(3);
      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);bestStrategy.push_back(2);
      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);bestStrategy.push_back(3);
      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);bestStrategy.push_back(3);
      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);bestStrategy.push_back(3);
      //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); bestStrategy.push_back(3);
      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); bestStrategy.push_back(3);
      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);bestStrategy.push_back(3);
      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);bestStrategy.push_back(3);
      //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);bestStrategy.push_back(3);
      //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);bestStrategy.push_back(3);
      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);bestStrategy.push_back(3);
      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);bestStrategy.push_back(3);
      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);bestStrategy.push_back(3);
      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);bestStrategy.push_back(3);
    } else {
      // Just testing one
      mpsName.push_back(empty.problemName());min.push_back(true);nRows.push_back(-1);
      nCols.push_back(-1);objValueTol.push_back(1.e-10);
      objValue.push_back(0.0);bestStrategy.push_back(0);
      int iTest;
      std::string alg;
      for (iTest=0;iTest<NUMBER_ALGORITHMS;iTest++) {
        ClpSolve solveOptions=setupForSolve(iTest,alg,0);
        printf("%d %s ",iTest,alg.c_str());
        if (switchOff[iTest]) 
          printf("skipped by user\n");
        else if(solveOptions.getSolveType()==ClpSolve::notImplemented)
          printf("skipped as not available\n");
        else
          printf("will be tested\n");
      }
    }

    double timeTaken =0.0;
    if( !barrierAvailable)
      switchOff[0]=1;
  // Loop once for each Mps File
    for (m=0; m<mpsName.size(); m++ ) {
      std::cerr <<"  processing mps file: " <<mpsName[m] 
		<<" (" <<m+1 <<" out of " <<mpsName.size() <<")" <<std::endl;

      ClpSimplex solutionBase=empty;
      std::string fn = dirNetlib+mpsName[m];
      if (!empty.numberRows()||algorithm<6) {
        // Read data mps file,
        CoinMpsIO mps;
	int nerrors=mps.readMps(fn.c_str(),"mps");
	if (nerrors) {
	  std::cerr << "Error " << nerrors << " when reading model from "
		    << fn.c_str() << "! Aborting tests.\n";
	  return 1;
	}
	solutionBase.loadProblem(*mps.getMatrixByCol(),mps.getColLower(),
                                 mps.getColUpper(),
                                 mps.getObjCoefficients(),
                                 mps.getRowLower(),mps.getRowUpper());
        
        solutionBase.setDblParam(ClpObjOffset,mps.objectiveOffset());
      } 
      
      // Runs through strategies
      if (algorithm==6||algorithm==7) {
        // algorithms tested are at top of file
        double testTime[NUMBER_ALGORITHMS];
        std::string alg[NUMBER_ALGORITHMS];
        int iTest;
        for (iTest=0;iTest<NUMBER_ALGORITHMS;iTest++) {
          ClpSolve solveOptions=setupForSolve(iTest,alg[iTest],1);
          if (solveOptions.getSolveType()!=ClpSolve::notImplemented) {
            double time1 = CoinCpuTime();
            ClpSimplex solution=solutionBase;
            if (solution.maximumSeconds()<0.0)
              solution.setMaximumSeconds(120.0);
	    if (doVector) {
	      ClpMatrixBase * matrix = solution.clpMatrix();
	      if (dynamic_cast< ClpPackedMatrix*>(matrix)) {
		ClpPackedMatrix * clpMatrix = dynamic_cast< ClpPackedMatrix*>(matrix);
		clpMatrix->makeSpecialColumnCopy();
	      }
	    }
            solution.initialSolve(solveOptions);
            double time2 = CoinCpuTime()-time1;
            testTime[iTest]=time2;
            printf("Took %g seconds - status %d\n",time2,solution.problemStatus());
            if (solution.problemStatus()) 
              testTime[iTest]=1.0e20;
          } else {
            testTime[iTest]=1.0e30;
          }
        }
        int iBest=-1;
        double dBest=1.0e10;
        printf("%s",fn.c_str());
        for (iTest=0;iTest<NUMBER_ALGORITHMS;iTest++) {
          if (testTime[iTest]<1.0e30) {
            printf(" %s %g",
                   alg[iTest].c_str(),testTime[iTest]);
            if (testTime[iTest]<dBest) {
              dBest=testTime[iTest];
              iBest=iTest;
            }
          }
        }
        printf("\n");
        if (iBest>=0)
          printf("Best strategy for %s is %s (%d) which takes %g seconds\n",
                 fn.c_str(),alg[iBest].c_str(),iBest,testTime[iBest]);
        else
          printf("No strategy finished in time limit\n");
        continue;
      }
      double time1 = CoinCpuTime();
      ClpSimplex solution=solutionBase;
#if 0
      solution.setOptimizationDirection(-1);
      {
	int j;
	double * obj = solution.objective();
	int n=solution.numberColumns();
	for (j=0;j<n;j++) 
	  obj[j] *= -1.0;
      }
#endif
      ClpSolve::SolveType method;
      ClpSolve::PresolveType presolveType;
      ClpSolve solveOptions;
      if (doPresolve)
        presolveType=ClpSolve::presolveOn;
      else
        presolveType=ClpSolve::presolveOff;
      solveOptions.setPresolveType(presolveType,5);
      std::string nameAlgorithm;
      if (algorithm!=5) {
        if (algorithm==0) {
          method=ClpSolve::useDual;
          nameAlgorithm="dual";
        } else if (algorithm==1) {
          method=ClpSolve::usePrimalorSprint;
          nameAlgorithm="primal";
        } else if (algorithm==3) {
          method=ClpSolve::automatic;
          nameAlgorithm="either";
        } else {
          nameAlgorithm="barrier-slow";
#ifdef UFL_BARRIER
          solveOptions.setSpecialOption(4,4);
          nameAlgorithm="barrier-UFL";
#endif
#ifdef WSSMP_BARRIER
          solveOptions.setSpecialOption(4,2);
          nameAlgorithm="barrier-WSSMP";
#endif
          method = ClpSolve::useBarrier;
        }
        solveOptions.setSolveType(method);
      } else {
        int iAlg = bestStrategy[m];
        int presolveOff=iAlg/100;
        iAlg=iAlg % 100;
        if( !barrierAvailable&&iAlg==0) {
          if (nRows[m]!=2172)
            iAlg = 5; // try primal
          else
            iAlg=3; // d2q06c
	}
        solveOptions=setupForSolve(iAlg,nameAlgorithm,0);
        if (presolveOff)
          solveOptions.setPresolveType(ClpSolve::presolveOff);
      }
      if (doVector) {
	ClpMatrixBase * matrix = solution.clpMatrix();
	if (dynamic_cast< ClpPackedMatrix*>(matrix)) {
	  ClpPackedMatrix * clpMatrix = dynamic_cast< ClpPackedMatrix*>(matrix);
	  clpMatrix->makeSpecialColumnCopy();
	}
      }
      solution.initialSolve(solveOptions);
      double time2 = CoinCpuTime()-time1;
      timeTaken += time2;
      printf("%s took %g seconds using algorithm %s\n",fn.c_str(),time2,nameAlgorithm.c_str());
      // Test objective solution value
      {
        double soln = solution.objectiveValue();
        CoinRelFltEq eq(objValueTol[m]);
        std::cerr <<soln <<",  " <<objValue[m] <<" diff "<<
	  soln-objValue[m]<<std::endl;
        if(!eq(soln,objValue[m]))
	  printf("** difference fails\n");
      }
    }
    printf("Total time %g seconds\n",timeTaken);
  }
  else {
    testingMessage( "***Skipped Testing on netlib    ***\n" );
    testingMessage( "***use -netlib to test class***\n" );
  }
  
  testingMessage( "All tests completed successfully\n" );
  return 0;
}

 
// Display message on stdout and stderr
void testingMessage( const char * const msg )
{
  std::cerr <<msg;
  //cout <<endl <<"*****************************************"
  //     <<endl <<msg <<endl;
}

//--------------------------------------------------------------------------
// test factorization methods and simplex method and simple barrier
void
ClpSimplexUnitTest(const std::string & dirSample)
{
  
  CoinRelFltEq eq(0.000001);

  {
    ClpSimplex solution;
  
    // matrix data
    //deliberate hiccup of 2 between 0 and 1
    CoinBigIndex start[5]={0,4,7,8,9};
    int length[5]={2,3,1,1,1};
    int rows[11]={0,2,-1,-1,0,1,2,0,1,2};
    double elements[11]={7.0,2.0,1.0e10,1.0e10,-2.0,1.0,-2.0,1,1,1};
    CoinPackedMatrix matrix(true,3,5,8,elements,rows,start,length);
    
    // rim data
    double objective[7]={-4.0,1.0,0.0,0.0,0.0,0.0,0.0};
    double rowLower[5]={14.0,3.0,3.0,1.0e10,1.0e10};
    double rowUpper[5]={14.0,3.0,3.0,-1.0e10,-1.0e10};
    double colLower[7]={0.0,0.0,0.0,0.0,0.0,0.0,0.0};
    double colUpper[7]={100.0,100.0,100.0,100.0,100.0,100.0,100.0};
    
    // basis 1
    int rowBasis1[3]={-1,-1,-1};
    int colBasis1[5]={1,1,-1,-1,1};
    solution.loadProblem(matrix,colLower,colUpper,objective,
			 rowLower,rowUpper);
    int i;
    solution.createStatus();
    for (i=0;i<3;i++) {
      if (rowBasis1[i]<0) {
	solution.setRowStatus(i,ClpSimplex::atLowerBound);
      } else {
	solution.setRowStatus(i,ClpSimplex::basic);
      }
    }
    for (i=0;i<5;i++) {
      if (colBasis1[i]<0) {
	solution.setColumnStatus(i,ClpSimplex::atLowerBound);
      } else {
	solution.setColumnStatus(i,ClpSimplex::basic);
      }
    }
    solution.setLogLevel(3+4+8+16+32);
    solution.primal();
    for (i=0;i<3;i++) {
      if (rowBasis1[i]<0) {
	solution.setRowStatus(i,ClpSimplex::atLowerBound);
      } else {
	solution.setRowStatus(i,ClpSimplex::basic);
      }
    }
    for (i=0;i<5;i++) {
      if (colBasis1[i]<0) {
	solution.setColumnStatus(i,ClpSimplex::atLowerBound);
      } else {
	solution.setColumnStatus(i,ClpSimplex::basic);
      }
    }
    // intricate stuff does not work with scaling
    solution.scaling(0);
    int returnCode = solution.factorize ( );
    assert(!returnCode);
    const double * colsol = solution.primalColumnSolution();
    const double * rowsol = solution.primalRowSolution();
    solution.getSolution(rowsol,colsol);
    double colsol1[5]={20.0/7.0,3.0,0.0,0.0,23.0/7.0};
    for (i=0;i<5;i++) {
      assert(eq(colsol[i],colsol1[i]));
    }
    // now feed in again without actually doing factorization
    ClpFactorization factorization2 = *solution.factorization();
    ClpSimplex solution2 = solution;
    solution2.setFactorization(factorization2);
    solution2.createStatus();
    for (i=0;i<3;i++) {
      if (rowBasis1[i]<0) {
	solution2.setRowStatus(i,ClpSimplex::atLowerBound);
      } else {
	solution2.setRowStatus(i,ClpSimplex::basic);
      }
    }
    for (i=0;i<5;i++) {
      if (colBasis1[i]<0) {
	solution2.setColumnStatus(i,ClpSimplex::atLowerBound);
      } else {
	solution2.setColumnStatus(i,ClpSimplex::basic);
      }
    }
    // intricate stuff does not work with scaling
    solution2.scaling(0);
    solution2.getSolution(rowsol,colsol);
    colsol = solution2.primalColumnSolution();
    rowsol = solution2.primalRowSolution();
    for (i=0;i<5;i++) {
      assert(eq(colsol[i],colsol1[i]));
    }
    solution2.setDualBound(0.1);
    solution2.dual();
    objective[2]=-1.0;
    objective[3]=-0.5;
    objective[4]=10.0;
    solution.dual();
    for (i=0;i<3;i++) {
      rowLower[i]=-1.0e20;
      colUpper[i+2]=0.0;
    }
    solution.setLogLevel(3);
    solution.dual();
    double rowObjective[]={1.0,0.5,-10.0};
    solution.loadProblem(matrix,colLower,colUpper,objective,
			 rowLower,rowUpper,rowObjective);
    solution.dual();
    solution.loadProblem(matrix,colLower,colUpper,objective,
			 rowLower,rowUpper,rowObjective);
    solution.primal();
  }
#ifndef COIN_NO_CLP_MESSAGE
  {    
    CoinMpsIO m;
    std::string fn = dirSample+"exmip1";
    m.readMps(fn.c_str(),"mps");
    ClpSimplex solution;
    solution.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
			 m.getObjCoefficients(),
			 m.getRowLower(),m.getRowUpper());
    solution.dual();
    // Test event handling
    MyEventHandler handler;
    solution.passInEventHandler(&handler);
    int numberRows=solution.numberRows();
    // make sure values pass has something to do
    for (int i=0;i<numberRows;i++)
      solution.setRowStatus(i,ClpSimplex::basic);
    solution.primal(1);
    assert (solution.secondaryStatus()==102); // Came out at end of pass
  }
  // Test Message handler
  {    
    CoinMpsIO m;
    std::string fn = dirSample+"exmip1";
    //fn = "Test/subGams4";
    m.readMps(fn.c_str(),"mps");
    ClpSimplex model;
    model.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
			 m.getObjCoefficients(),
			 m.getRowLower(),m.getRowUpper());
    // Message handler
    MyMessageHandler messageHandler(&model);
    std::cout<<"Testing derived message handler"<<std::endl;
    model.passInMessageHandler(&messageHandler);
    model.messagesPointer()->setDetailMessage(1,102);
    model.setFactorizationFrequency(10);
    model.primal();
    model.primal(0,3);
    model.setObjCoeff(3,-2.9473684210526314);
    model.primal(0,3);
    // Write saved solutions
    int nc = model.getNumCols();
    int s; 
    std::deque<StdVectorDouble> fep = messageHandler.getFeasibleExtremePoints();
    int numSavedSolutions = fep.size();
    for ( s=0; s<numSavedSolutions; ++s ) {
      const StdVectorDouble & solnVec = fep[s];
      for ( int c=0; c<nc; ++c ) {
        if (fabs(solnVec[c])>1.0e-8)
          std::cout <<"Saved Solution: " <<s <<" ColNum: " <<c <<" Value: " <<solnVec[c] <<std::endl;
      }
    }
    // Solve again without scaling
    // and maximize then minimize
    messageHandler.clearFeasibleExtremePoints();
    model.scaling(0);
    model.setOptimizationDirection(-1);
    model.primal();
    model.setOptimizationDirection(1);
    model.primal();
    fep = messageHandler.getFeasibleExtremePoints();
    numSavedSolutions = fep.size();
    for ( s=0; s<numSavedSolutions; ++s ) {
      const StdVectorDouble & solnVec = fep[s];
      for ( int c=0; c<nc; ++c ) {
        if (fabs(solnVec[c])>1.0e-8)
          std::cout <<"Saved Solution: " <<s <<" ColNum: " <<c <<" Value: " <<solnVec[c] <<std::endl;
      }
    }
  }
#endif
  // Test dual ranging
  {    
    CoinMpsIO m;
    std::string fn = dirSample+"exmip1";
    m.readMps(fn.c_str(),"mps");
    ClpSimplex model;
    model.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
			 m.getObjCoefficients(),
			 m.getRowLower(),m.getRowUpper());
    model.primal();
    int which[13] = {0,1,2,3,4,5,6,7,8,9,10,11,12};
    double costIncrease[13];
    int sequenceIncrease[13];
    double costDecrease[13];
    int sequenceDecrease[13];
    // ranging
    model.dualRanging(13,which,costIncrease,sequenceIncrease,
		      costDecrease,sequenceDecrease);
    int i;
    for ( i=0;i<13;i++)
      printf("%d increase %g %d, decrease %g %d\n",
	     i,costIncrease[i],sequenceIncrease[i],
	     costDecrease[i],sequenceDecrease[i]);
    assert (fabs(costDecrease[3])<1.0e-4);
    assert (fabs(costIncrease[7]-1.0)<1.0e-4);
    model.setOptimizationDirection(-1);
    {
      int j;
      double * obj = model.objective();
      int n=model.numberColumns();
      for (j=0;j<n;j++) 
	obj[j] *= -1.0;
    }
    double costIncrease2[13];
    int sequenceIncrease2[13];
    double costDecrease2[13];
    int sequenceDecrease2[13];
    // ranging
    model.dualRanging(13,which,costIncrease2,sequenceIncrease2,
		      costDecrease2,sequenceDecrease2);
    for (i=0;i<13;i++) {
      assert (fabs(costIncrease[i]-costDecrease2[i])<1.0e-6);
      assert (fabs(costDecrease[i]-costIncrease2[i])<1.0e-6);
      assert (sequenceIncrease[i]==sequenceDecrease2[i]);
      assert (sequenceDecrease[i]==sequenceIncrease2[i]);
    }
    // Now delete all rows and see what happens
    model.deleteRows(model.numberRows(),which);
    model.primal();
    // ranging
    if (!model.dualRanging(8,which,costIncrease,sequenceIncrease,
                           costDecrease,sequenceDecrease)) {
      for (i=0;i<8;i++) {
        printf("%d increase %g %d, decrease %g %d\n",
               i,costIncrease[i],sequenceIncrease[i],
               costDecrease[i],sequenceDecrease[i]);
      }
    }
  }
  // Test primal ranging
  {    
    CoinMpsIO m;
    std::string fn = dirSample+"exmip1";
    m.readMps(fn.c_str(),"mps");
    ClpSimplex model;
    model.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
			 m.getObjCoefficients(),
			 m.getRowLower(),m.getRowUpper());
    model.primal();
    int which[13] = {0,1,2,3,4,5,6,7,8,9,10,11,12};
    double valueIncrease[13];
    int sequenceIncrease[13];
    double valueDecrease[13];
    int sequenceDecrease[13];
    // ranging
    model.primalRanging(13,which,valueIncrease,sequenceIncrease,
		      valueDecrease,sequenceDecrease);
    int i;
    for ( i=0;i<13;i++)
      printf("%d increase %g %d, decrease %g %d\n",
	     i,valueIncrease[i],sequenceIncrease[i],
	     valueDecrease[i],sequenceDecrease[i]);
    assert (fabs(valueIncrease[3]-0.642857)<1.0e-4);
    assert (fabs(valueIncrease[8]-2.95113)<1.0e-4);
#if 0
    // out until I find optimization bug
    // Test parametrics
    ClpSimplexOther * model2 = (ClpSimplexOther *) (&model);
    double rhs[]={ 1.0,2.0,3.0,4.0,5.0};
    double endingTheta=1.0;
    model2->scaling(0);
    model2->setLogLevel(63);
    model2->parametrics(0.0,endingTheta,0.1,
                        NULL,NULL,rhs,rhs,NULL);
#endif
  }
  // Test binv etc
  {    
    /* 
       Wolsey : Page 130
       max 4x1 -  x2
       7x1 - 2x2    <= 14
       x2    <= 3
       2x1 - 2x2    <= 3
       x1 in Z+, x2 >= 0

       note slacks are -1 in Clp so signs may be different
    */
    
    int n_cols = 2;
    int n_rows = 3;
    
    double obj[2] = {-4.0, 1.0};
    double collb[2] = {0.0, 0.0};
    double colub[2] = {COIN_DBL_MAX, COIN_DBL_MAX};
    double rowlb[3] = {-COIN_DBL_MAX, -COIN_DBL_MAX, -COIN_DBL_MAX};
    double rowub[3] = {14.0, 3.0, 3.0};
    
    int rowIndices[5] =  {0,     2,    0,    1,    2};
    int colIndices[5] =  {0,     0,    1,    1,    1};
    double elements[5] = {7.0, 2.0, -2.0,  1.0, -2.0};
    CoinPackedMatrix M(true, rowIndices, colIndices, elements, 5);

    ClpSimplex model;
    model.loadProblem(M, collb, colub, obj, rowlb, rowub);
    model.dual(0,1); // keep factorization 
    
    //check that the tableau matches wolsey (B-1 A)
    // slacks in second part of binvA
    double * binvA = (double*) malloc((n_cols+n_rows) * sizeof(double));
    
    printf("B-1 A by row\n");
    int i;
    for( i = 0; i < n_rows; i++){
      model.getBInvARow(i, binvA,binvA+n_cols);
      printf("row: %d -> ",i);
      for(int j=0; j < n_cols+n_rows; j++){
	printf("%g, ", binvA[j]);
      }
      printf("\n");
    }
    // See if can re-use factorization AND arrays
    model.primal(0,3+4); // keep factorization
    // And do by column
    printf("B-1 A by column\n");
    for( i = 0; i < n_rows+n_cols; i++){
      model.getBInvACol(i, binvA);
      printf("column: %d -> ",i);
      for(int j=0; j < n_rows; j++){
	printf("%g, ", binvA[j]);
      }
      printf("\n");
    }
    /* Do twice -
       without and with scaling
    */
    // set scaling off
    model.scaling(0);
    for (int iPass=0;iPass<2;iPass++) {
      model.primal(0,3+4); // keep factorization
      const double * rowScale = model.rowScale();
      const double * columnScale = model.columnScale();
      if (!iPass)
        assert (!rowScale);
      else
        assert (rowScale); // only true for this example
      /* has to be exactly correct as in OsiClpsolverInterface.cpp
         (also redo each pass as may change
      */
      printf("B-1 A");
      for( i = 0; i < n_rows; i++){
        model.getBInvARow(i, binvA,binvA+n_cols);
        printf("\nrow: %d -> ",i);
        int j;
        // First columns
        for(j=0; j < n_cols; j++){
	  if (binvA[j]) {
	    printf("(%d %g), ", j, binvA[j]);
	  }
        }
        // now rows
        for(j=0; j < n_rows; j++){
	  if (binvA[j+n_cols]) {
	    printf("(%d %g), ", j+n_cols, binvA[j+n_cols]);
	  }
        }
      }
      printf("\n");
      printf("And by column (trickier)");
      const int * pivotVariable = model.pivotVariable();
      for( i = 0; i < n_cols+n_rows; i++){
        model.getBInvACol(i, binvA);
        printf("\ncolumn: %d -> ",i);
        for(int j=0; j < n_rows; j++){
	  if (binvA[j]) {
	    // need to know pivot variable for +1/-1 (slack) and row/column scaling
	    int pivot = pivotVariable[j];
	    if (pivot<n_cols) {
	      // scaled coding is in just in case 
	      if (!columnScale) {
		printf("(%d %g), ", j, binvA[j]);
	      } else {
		printf("(%d %g), ", j, binvA[j]*columnScale[pivot]);
	      }
	    } else {
	      if (!rowScale) {
		printf("(%d %g), ", j, binvA[j]);
	      } else {
		printf("(%d %g), ", j, binvA[j]/rowScale[pivot-n_cols]);
	      }
	    }
	  }
        }
      }
      printf("\n");
      printf("binvrow");
      for( i = 0; i < n_rows; i++){
	model.getBInvRow(i, binvA);
	printf("\nrow: %d -> ",i);
	int j;
	for (j=0;j<n_rows;j++) {
	  if (binvA[j])
	    printf("(%d %g), ", j, binvA[j]);
	}
      }
      printf("\n");
      printf("And by column ");
      for( i = 0; i < n_rows; i++){
	model.getBInvCol(i, binvA);
	printf("\ncol: %d -> ",i);
	int j;
	for (j=0;j<n_rows;j++) {
	  if (binvA[j])
	    printf("(%d %g), ", j, binvA[j]);
	}
      }
      printf("\n");
      // now deal with next pass
      if (!iPass) {
	// get scaling for testing
	model.scaling(1);
      }
    }
    free(binvA);
    model.setColUpper(1,2.0);
    model.dual(0,2+4); // use factorization and arrays
    model.dual(0,2); // hopefully will not use factorization
    model.primal(0,3+4); // keep factorization
    // but say basis has changed
    model.setWhatsChanged(model.whatsChanged()&(~512));
    model.dual(0,2); // hopefully will not use factorization
  }
  // test steepest edge
  {    
    CoinMpsIO m;
    std::string fn = dirSample+"finnis";
    int returnCode = m.readMps(fn.c_str(),"mps");
    if (returnCode) {
      // probable cause is that gz not there
      fprintf(stderr,"Unable to open finnis.mps in %s\n", dirSample.c_str());
      fprintf(stderr,"Most probable cause is finnis.mps is gzipped i.e. finnis.mps.gz and libz has not been activated\n");
      fprintf(stderr,"Either gunzip files or edit Makefiles/Makefile.location to get libz\n");
      exit(999);
    }
    ClpModel model;
    model.loadProblem(*m.getMatrixByCol(),m.getColLower(),
		    m.getColUpper(),
		    m.getObjCoefficients(),
		    m.getRowLower(),m.getRowUpper());
    ClpSimplex solution(model);

    solution.scaling(1); 
    solution.setDualBound(1.0e8);
    //solution.factorization()->maximumPivots(1);
    //solution.setLogLevel(3);
    solution.setDualTolerance(1.0e-7);
    // set objective sense,
    ClpDualRowSteepest steep;
    solution.setDualRowPivotAlgorithm(steep);
    solution.setDblParam(ClpObjOffset,m.objectiveOffset());
    solution.dual();
  }
  // test normal solution
  {    
    CoinMpsIO m;
    std::string fn = dirSample+"afiro";
    m.readMps(fn.c_str(),"mps");
    ClpSimplex solution;
    ClpModel model;
    // do twice - without and with scaling
    int iPass;
    for (iPass=0;iPass<2;iPass++) {
      // explicit row objective for testing
      int nr = m.getNumRows();
      double * rowObj = new double[nr];
      CoinFillN(rowObj,nr,0.0);
      model.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
		      m.getObjCoefficients(),
		      m.getRowLower(),m.getRowUpper(),rowObj);
      delete [] rowObj;
      solution = ClpSimplex(model);
      if (iPass) {
	solution.scaling();
      }
      solution.dual();
      solution.dual();
      // test optimal
      assert (solution.status()==0);
      int numberColumns = solution.numberColumns();
      int numberRows = solution.numberRows();
      CoinPackedVector colsol(numberColumns,solution.primalColumnSolution());
      double * objective = solution.objective();
      double objValue = colsol.dotProduct(objective);
      CoinRelFltEq eq(1.0e-8);
      assert(eq(objValue,-4.6475314286e+02));
      // Test auxiliary model
      //solution.scaling(0);
      solution.auxiliaryModel(63-2); // bounds may change
      solution.dual();
      solution.primal();
      solution.allSlackBasis();
      solution.dual();
      assert(eq(solution.objectiveValue(),-4.6475314286e+02));
      solution.auxiliaryModel(-1);
      solution.dual();
      assert(eq(solution.objectiveValue(),-4.6475314286e+02));
      double * lower = solution.columnLower();
      double * upper = solution.columnUpper();
      double * sol = solution.primalColumnSolution();
      double * result = new double[numberColumns];
      CoinFillN ( result, numberColumns,0.0);
      solution.matrix()->transposeTimes(solution.dualRowSolution(), result);
      int iRow , iColumn;
      // see if feasible and dual feasible
      for (iColumn=0;iColumn<numberColumns;iColumn++) {
	double value = sol[iColumn];
	assert(value<upper[iColumn]+1.0e-8);
	assert(value>lower[iColumn]-1.0e-8);
	value = objective[iColumn]-result[iColumn];
	assert (value>-1.0e-5);
	if (sol[iColumn]>1.0e-5)
	  assert (value<1.0e-5);
      }
      delete [] result;
      result = new double[numberRows];
      CoinFillN ( result, numberRows,0.0);
      solution.matrix()->times(colsol, result);
      lower = solution.rowLower();
      upper = solution.rowUpper();
      sol = solution.primalRowSolution();
      for (iRow=0;iRow<numberRows;iRow++) {
	double value = result[iRow];
	assert(eq(value,sol[iRow]));
	assert(value<upper[iRow]+1.0e-8);
	assert(value>lower[iRow]-1.0e-8);
      }
      delete [] result;
      // test row objective
      double * rowObjective = solution.rowObjective();
      CoinDisjointCopyN(solution.dualRowSolution(),numberRows,rowObjective);
      CoinDisjointCopyN(solution.dualColumnSolution(),numberColumns,objective);
      // this sets up all slack basis
      solution.createStatus();
      solution.dual();
      CoinFillN(rowObjective,numberRows,0.0);
      CoinDisjointCopyN(m.getObjCoefficients(),numberColumns,objective);
      solution.dual();
    }
  }
  // test unbounded
  {    
    CoinMpsIO m;
    std::string fn = dirSample+"brandy";
    m.readMps(fn.c_str(),"mps");
    ClpSimplex solution;
    // do twice - without and with scaling
    int iPass;
    for (iPass=0;iPass<2;iPass++) {
      solution.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
		      m.getObjCoefficients(),
		      m.getRowLower(),m.getRowUpper());
      if (iPass)
	solution.scaling();
      solution.setOptimizationDirection(-1);
      // test unbounded and ray
#ifdef DUAL
      solution.setDualBound(100.0);
      solution.dual();
#else
      solution.primal();
#endif
      assert (solution.status()==2);
      int numberColumns = solution.numberColumns();
      int numberRows = solution.numberRows();
      double * lower = solution.columnLower();
      double * upper = solution.columnUpper();
      double * sol = solution.primalColumnSolution();
      double * ray = solution.unboundedRay();
      double * objective = solution.objective();
      double objChange=0.0;
      int iRow , iColumn;
      // make sure feasible and columns form ray
      for (iColumn=0;iColumn<numberColumns;iColumn++) {
	double value = sol[iColumn];
	assert(value<upper[iColumn]+1.0e-8);
	assert(value>lower[iColumn]-1.0e-8);
	value = ray[iColumn];
	if (value>0.0)
	  assert(upper[iColumn]>1.0e30);
	else if (value<0.0)
	  assert(lower[iColumn]<-1.0e30);
	objChange += value*objective[iColumn];
      }
      // make sure increasing objective
      assert(objChange>0.0);
      double * result = new double[numberRows];
      CoinFillN ( result, numberRows,0.0);
      solution.matrix()->times(sol, result);
      lower = solution.rowLower();
      upper = solution.rowUpper();
      sol = solution.primalRowSolution();
      for (iRow=0;iRow<numberRows;iRow++) {
	double value = result[iRow];
	assert(eq(value,sol[iRow]));
	assert(value<upper[iRow]+2.0e-8);
	assert(value>lower[iRow]-2.0e-8);
      }
      CoinFillN ( result, numberRows,0.0);
      solution.matrix()->times(ray, result);
      // there may be small differences (especially if scaled)
      for (iRow=0;iRow<numberRows;iRow++) {
	double value = result[iRow];
	if (value>1.0e-8)
	  assert(upper[iRow]>1.0e30);
	else if (value<-1.0e-8)
	  assert(lower[iRow]<-1.0e30);
      }
      delete [] result;
      delete [] ray;
    }
  }
  // test infeasible
  {    
    CoinMpsIO m;
    std::string fn = dirSample+"brandy";
    m.readMps(fn.c_str(),"mps");
    ClpSimplex solution;
    // do twice - without and with scaling
    int iPass;
    for (iPass=0;iPass<2;iPass++) {
      solution.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
		      m.getObjCoefficients(),
		      m.getRowLower(),m.getRowUpper());
      if (iPass)
	solution.scaling();
      // test infeasible and ray
      solution.columnUpper()[0]=0.0;
#ifdef DUAL
      solution.setDualBound(100.0);
      solution.dual();
#else
      solution.primal();
#endif
      assert (solution.status()==1);
      int numberColumns = solution.numberColumns();
      int numberRows = solution.numberRows();
      double * lower = solution.rowLower();
      double * upper = solution.rowUpper();
      double * ray = solution.infeasibilityRay();
      assert(ray);
      // construct proof of infeasibility
      int iRow , iColumn;
      double lo=0.0,up=0.0;
      int nl=0,nu=0;
      for (iRow=0;iRow<numberRows;iRow++) {
	if (lower[iRow]>-1.0e20) {
	  lo += ray[iRow]*lower[iRow];
	} else {
	  if (ray[iRow]>1.0e-8) 
	    nl++;
	}
	if (upper[iRow]<1.0e20) {
	  up += ray[iRow]*upper[iRow];
	} else {
	  if (ray[iRow]>1.0e-8) 
	    nu++;
	}
      }
      if (nl)
	lo=-1.0e100;
      if (nu)
	up=1.0e100;
      double * result = new double[numberColumns];
      double lo2=0.0,up2=0.0;
      CoinFillN ( result, numberColumns,0.0);
      solution.matrix()->transposeTimes(ray, result);
      lower = solution.columnLower();
      upper = solution.columnUpper();
      nl=nu=0;
      for (iColumn=0;iColumn<numberColumns;iColumn++) {
	if (result[iColumn]>1.0e-8) {
	  if (lower[iColumn]>-1.0e20)
	    lo2 += result[iColumn]*lower[iColumn];
	  else
	    nl++;
	  if (upper[iColumn]<1.0e20)
	    up2 += result[iColumn]*upper[iColumn];
	  else
	    nu++;
	} else if (result[iColumn]<-1.0e-8) {
	  if (lower[iColumn]>-1.0e20)
	    up2 += result[iColumn]*lower[iColumn];
	  else
	    nu++;
	  if (upper[iColumn]<1.0e20)
	    lo2 += result[iColumn]*upper[iColumn];
	  else
	    nl++;
	}
      }
      if (nl)
	lo2=-1.0e100;
      if (nu)
	up2=1.0e100;
      // make sure inconsistency
      assert(lo2>up||up2<lo);
      delete [] result;
      delete [] ray;
    }
  }
  // test delete and add
  {    
    CoinMpsIO m;
    std::string fn = dirSample+"brandy";
    m.readMps(fn.c_str(),"mps");
    ClpSimplex solution;
    solution.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
			 m.getObjCoefficients(),
			 m.getRowLower(),m.getRowUpper());
    solution.dual();
    CoinRelFltEq eq(1.0e-8);
    assert(eq(solution.objectiveValue(),1.5185098965e+03));

    int numberColumns = solution.numberColumns();
    int numberRows = solution.numberRows();
    double * saveObj = new double [numberColumns];
    double * saveLower = new double[numberRows+numberColumns];
    double * saveUpper = new double[numberRows+numberColumns];
    int * which = new int [numberRows+numberColumns];

    int numberElements = m.getMatrixByCol()->getNumElements();
    int * starts = new int[numberRows+numberColumns];
    int * index = new int[numberElements];
    double * element = new double[numberElements];

    const CoinBigIndex * startM;
    const int * lengthM;
    const int * indexM;
    const double * elementM;

    int n,nel;

    // delete non basic columns
    n=0;
    nel=0;
    int iRow , iColumn;
    const double * lower = m.getColLower();
    const double * upper = m.getColUpper();
    const double * objective = m.getObjCoefficients();
    startM = m.getMatrixByCol()->getVectorStarts();
    lengthM = m.getMatrixByCol()->getVectorLengths();
    indexM = m.getMatrixByCol()->getIndices();
    elementM = m.getMatrixByCol()->getElements();
    starts[0]=0;
    for (iColumn=0;iColumn<numberColumns;iColumn++) {
      if (solution.getColumnStatus(iColumn)!=ClpSimplex::basic) {
	saveObj[n]=objective[iColumn];
	saveLower[n]=lower[iColumn];
	saveUpper[n]=upper[iColumn];
	int j;
	for (j=startM[iColumn];j<startM[iColumn]+lengthM[iColumn];j++) {
	  index[nel]=indexM[j];
	  element[nel++]=elementM[j];
	}
	which[n++]=iColumn;
	starts[n]=nel;
      }
    }
    solution.deleteColumns(n,which);
    solution.dual();
    // Put back
    solution.addColumns(n,saveLower,saveUpper,saveObj,
			starts,index,element);
    solution.dual();
    assert(eq(solution.objectiveValue(),1.5185098965e+03));
    // Delete all columns and add back
    n=0;
    nel=0;
    starts[0]=0;
    lower = m.getColLower();
    upper = m.getColUpper();
    objective = m.getObjCoefficients();
    for (iColumn=0;iColumn<numberColumns;iColumn++) {
      saveObj[n]=objective[iColumn];
      saveLower[n]=lower[iColumn];
      saveUpper[n]=upper[iColumn];
      int j;
      for (j=startM[iColumn];j<startM[iColumn]+lengthM[iColumn];j++) {
	index[nel]=indexM[j];
	element[nel++]=elementM[j];
      }
      which[n++]=iColumn;
      starts[n]=nel;
    }
    solution.deleteColumns(n,which);
    solution.dual();
    // Put back
    solution.addColumns(n,saveLower,saveUpper,saveObj,
			starts,index,element);
    solution.dual();
    assert(eq(solution.objectiveValue(),1.5185098965e+03));

    // reload with original
    solution.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
			 m.getObjCoefficients(),
			 m.getRowLower(),m.getRowUpper());
    // delete half rows
    n=0;
    nel=0;
    lower = m.getRowLower();
    upper = m.getRowUpper();
    startM = m.getMatrixByRow()->getVectorStarts();
    lengthM = m.getMatrixByRow()->getVectorLengths();
    indexM = m.getMatrixByRow()->getIndices();
    elementM = m.getMatrixByRow()->getElements();
    starts[0]=0;
    for (iRow=0;iRow<numberRows;iRow++) {
      if ((iRow&1)==0) {
	saveLower[n]=lower[iRow];
	saveUpper[n]=upper[iRow];
	int j;
	for (j=startM[iRow];j<startM[iRow]+lengthM[iRow];j++) {
	  index[nel]=indexM[j];
	  element[nel++]=elementM[j];
	}
	which[n++]=iRow;
	starts[n]=nel;
      }
    }
    solution.deleteRows(n,which);
    solution.dual();
    // Put back
    solution.addRows(n,saveLower,saveUpper,
			starts,index,element);
    solution.dual();
    assert(eq(solution.objectiveValue(),1.5185098965e+03));
    solution.writeMps("yy.mps");
    // Delete all rows
    n=0;
    nel=0;
    lower = m.getRowLower();
    upper = m.getRowUpper();
    starts[0]=0;
    for (iRow=0;iRow<numberRows;iRow++) {
      saveLower[n]=lower[iRow];
      saveUpper[n]=upper[iRow];
      int j;
      for (j=startM[iRow];j<startM[iRow]+lengthM[iRow];j++) {
	index[nel]=indexM[j];
	element[nel++]=elementM[j];
      }
      which[n++]=iRow;
      starts[n]=nel;
    }
    solution.deleteRows(n,which);
    solution.dual();
    // Put back
    solution.addRows(n,saveLower,saveUpper,
			starts,index,element);
    solution.dual();
    solution.writeMps("xx.mps");
    assert(eq(solution.objectiveValue(),1.5185098965e+03));
    // Zero out status array to give some interest
    memset(solution.statusArray()+numberColumns,0,numberRows);
    solution.primal(1);
    assert(eq(solution.objectiveValue(),1.5185098965e+03));
    // Delete all columns and rows
    n=0;
    for (iColumn=0;iColumn<numberColumns;iColumn++) {
      which[n++]=iColumn;
      starts[n]=nel;
    }
    solution.deleteColumns(n,which);
    n=0;
    for (iRow=0;iRow<numberRows;iRow++) {
      which[n++]=iRow;
      starts[n]=nel;
    }
    solution.deleteRows(n,which);

    delete [] saveObj;
    delete [] saveLower;
    delete [] saveUpper;
    delete [] which;
    delete [] starts;
    delete [] index;
    delete [] element;
  }
#if 1
  // Test barrier
  {
    CoinMpsIO m;
    std::string fn = dirSample+"exmip1";
    m.readMps(fn.c_str(),"mps");
    ClpInterior solution;
    solution.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
			 m.getObjCoefficients(),
			 m.getRowLower(),m.getRowUpper());
    solution.primalDual();
  }
#endif
  // test network 
#define QUADRATIC
  if (1) {    
    std::string fn = dirSample+"input.130";
    int numberColumns;
    int numberRows;
    
    FILE * fp = fopen(fn.c_str(),"r");
    if (!fp) {
      // Try in Data/Sample
      fn = "Data/Sample/input.130";
      fp = fopen(fn.c_str(),"r");
    }
    if (!fp) {
      fprintf(stderr,"Unable to open file input.130 in dirSample or Data/Sample directory\n");
    } else {
      int problem;
      char temp[100];
      // read and skip 
      fscanf(fp,"%s",temp);
      assert (!strcmp(temp,"BEGIN"));
      fscanf(fp,"%*s %*s %d %d %*s %*s %d %*s",&problem, &numberRows, 
	     &numberColumns);
      // scan down to SUPPLY
      while (fgets(temp,100,fp)) {
	if (!strncmp(temp,"SUPPLY",6))
	  break;
      }
      if (strncmp(temp,"SUPPLY",6)) {
	fprintf(stderr,"Unable to find SUPPLY\n");
	exit(2);
      }
      // get space for rhs
      double * lower = new double[numberRows];
      double * upper = new double[numberRows];
      int i;
      for (i=0;i<numberRows;i++) {
	lower[i]=0.0;
	upper[i]=0.0;
      }
      // ***** Remember to convert to C notation
      while (fgets(temp,100,fp)) {
	int row;
	int value;
	if (!strncmp(temp,"ARCS",4))
	  break;
	sscanf(temp,"%d %d",&row,&value);
	upper[row-1]=-value;
	lower[row-1]=-value;
      }
      if (strncmp(temp,"ARCS",4)) {
	fprintf(stderr,"Unable to find ARCS\n");
	exit(2);
      }
      // number of columns may be underestimate
      int * head = new int[2*numberColumns];
      int * tail = new int[2*numberColumns];
      int * ub = new int[2*numberColumns];
      int * cost = new int[2*numberColumns];
      // ***** Remember to convert to C notation
      numberColumns=0;
      while (fgets(temp,100,fp)) {
	int iHead;
	int iTail;
	int iUb;
	int iCost;
	if (!strncmp(temp,"DEMAND",6))
	  break;
	sscanf(temp,"%d %d %d %d",&iHead,&iTail,&iCost,&iUb);
	iHead--;
	iTail--;
	head[numberColumns]=iHead;
	tail[numberColumns]=iTail;
	ub[numberColumns]=iUb;
	cost[numberColumns]=iCost;
	numberColumns++;
      }
      if (strncmp(temp,"DEMAND",6)) {
	fprintf(stderr,"Unable to find DEMAND\n");
	exit(2);
      }
      // ***** Remember to convert to C notation
      while (fgets(temp,100,fp)) {
	int row;
	int value;
	if (!strncmp(temp,"END",3))
	  break;
	sscanf(temp,"%d %d",&row,&value);
	upper[row-1]=value;
	lower[row-1]=value;
      }
      if (strncmp(temp,"END",3)) {
	fprintf(stderr,"Unable to find END\n");
	exit(2);
      }
      printf("Problem %d has %d rows and %d columns\n",problem,
	     numberRows,numberColumns);
      fclose(fp);
      ClpSimplex  model;
      // now build model
      
      double * objective =new double[numberColumns];
      double * lowerColumn = new double[numberColumns];
      double * upperColumn = new double[numberColumns];
      
      double * element = new double [2*numberColumns];
      CoinBigIndex * start = new CoinBigIndex [numberColumns+1];
      int * row = new int[2*numberColumns];
      start[numberColumns]=2*numberColumns;
      for (i=0;i<numberColumns;i++) {
	start[i]=2*i;
	element[2*i]=-1.0;
	element[2*i+1]=1.0;
	row[2*i]=head[i];
	row[2*i+1]=tail[i];
	lowerColumn[i]=0.0;
	upperColumn[i]=ub[i];
	objective[i]=cost[i];
      }
      // Create Packed Matrix
      CoinPackedMatrix matrix;
      int * lengths = NULL;
      matrix.assignMatrix(true,numberRows,numberColumns,
			  2*numberColumns,element,row,start,lengths);
      // load model
      model.loadProblem(matrix,
			lowerColumn,upperColumn,objective,
			lower,upper);
      model.factorization()->maximumPivots(200+model.numberRows()/100);
      model.createStatus();
      double time1 = CoinCpuTime();
      model.dual();
      std::cout<<"Network problem, ClpPackedMatrix took "<<CoinCpuTime()-time1<<" seconds"<<std::endl;
      ClpPlusMinusOneMatrix * plusMinus = new ClpPlusMinusOneMatrix(matrix);
      assert (plusMinus->getIndices()); // would be zero if not +- one
      //ClpPlusMinusOneMatrix *plusminus_matrix;

      //plusminus_matrix = new ClpPlusMinusOneMatrix;

      //plusminus_matrix->passInCopy(numberRows, numberColumns, true, plusMinus->getMutableIndices(),
      //                         plusMinus->startPositive(),plusMinus->startNegative());
      model.loadProblem(*plusMinus,
			lowerColumn,upperColumn,objective,
			lower,upper);
      //model.replaceMatrix( plusminus_matrix , true);
      delete plusMinus;
      //model.createStatus();
      //model.initialSolve();
      //model.writeMps("xx.mps");
      
      model.factorization()->maximumPivots(200+model.numberRows()/100);
      model.createStatus();
      time1 = CoinCpuTime();
      model.dual();
      std::cout<<"Network problem, ClpPlusMinusOneMatrix took "<<CoinCpuTime()-time1<<" seconds"<<std::endl;
      ClpNetworkMatrix network(numberColumns,head,tail);
      model.loadProblem(network,
			lowerColumn,upperColumn,objective,
			lower,upper);
      
      model.factorization()->maximumPivots(200+model.numberRows()/100);
      model.createStatus();
      time1 = CoinCpuTime();
      model.dual();
      std::cout<<"Network problem, ClpNetworkMatrix took "<<CoinCpuTime()-time1<<" seconds"<<std::endl;
      delete [] lower;
      delete [] upper;
      delete [] head;
      delete [] tail;
      delete [] ub;
      delete [] cost;
      delete [] objective;
      delete [] lowerColumn;
      delete [] upperColumn;
    }
  }
#ifdef QUADRATIC
  // Test quadratic to solve linear
  if (1) {    
    CoinMpsIO m;
    std::string fn = dirSample+"exmip1";
    m.readMps(fn.c_str(),"mps");
    ClpSimplex solution;
    solution.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
			 m.getObjCoefficients(),
			 m.getRowLower(),m.getRowUpper());
    //solution.dual();
    // get quadratic part
    int numberColumns=solution.numberColumns();
    int * start=new int [numberColumns+1];
    int * column = new int[numberColumns];
    double * element = new double[numberColumns];
    int i;
    start[0]=0;
    int n=0;
    int kk=numberColumns-1;
    int kk2=numberColumns-1;
    for (i=0;i<numberColumns;i++) {
      if (i>=kk) {
	column[n]=i;
	if (i>=kk2)
	  element[n]=1.0e-1;
	else
	  element[n]=0.0;
	n++;
      }
      start[i+1]=n;
    }
    // Load up objective
    solution.loadQuadraticObjective(numberColumns,start,column,element);
    delete [] start;
    delete [] column;
    delete [] element;
    //solution.quadraticSLP(50,1.0e-4);
    CoinRelFltEq eq(1.0e-4);
    //assert(eq(objValue,820.0));
    //solution.setLogLevel(63);
    solution.primal();
    printSol(solution);
    //assert(eq(objValue,3.2368421));
    //exit(77);
  }
  // Test quadratic
  if (1) {    
    CoinMpsIO m;
    std::string fn = dirSample+"share2qp";
    //fn = "share2qpb";
    m.readMps(fn.c_str(),"mps");
    ClpSimplex model;
    model.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
			 m.getObjCoefficients(),
			 m.getRowLower(),m.getRowUpper());
    model.dual();
    // get quadratic part
    int * start=NULL;
    int * column = NULL;
    double * element = NULL;
    m.readQuadraticMps(NULL,start,column,element,2);
    int column2[200];
    double element2[200];
    int start2[80];
    int j;
    start2[0]=0;
    int nel=0;
    bool good=false;
    for (j=0;j<79;j++) {
      if (start[j]==start[j+1]) {
	column2[nel]=j;
	element2[nel]=0.0;
	nel++;
      } else {
	int i;
	for (i=start[j];i<start[j+1];i++) {
	  column2[nel]=column[i];
	  element2[nel++]=element[i];
	}
      }
      start2[j+1]=nel;
    }
    // Load up objective
    if (good)
      model.loadQuadraticObjective(model.numberColumns(),start2,column2,element2);
    else
      model.loadQuadraticObjective(model.numberColumns(),start,column,element);
    delete [] start;
    delete [] column;
    delete [] element;
    int numberColumns=model.numberColumns();
    model.scaling(0);
#if 0
    model.nonlinearSLP(50,1.0e-4);
#else
    // Get feasible
    ClpObjective * saveObjective = model.objectiveAsObject()->clone();
    ClpLinearObjective zeroObjective(NULL,numberColumns);
    model.setObjective(&zeroObjective);
    model.dual();
    model.setObjective(saveObjective);
    delete saveObjective;
#endif
    //model.setLogLevel(63);
    //exit(77);
    model.setFactorizationFrequency(10);
    model.primal();
    printSol(model);
    double objValue = model.getObjValue();
    CoinRelFltEq eq(1.0e-4);
    assert(eq(objValue,-400.92));
    // and again for barrier
    model.barrier(false);
    //printSol(model);
    model.allSlackBasis();
    model.primal();
    //printSol(model);
  }
  if (0) {    
    CoinMpsIO m;
    std::string fn = "./beale";
    //fn = "./jensen";
    m.readMps(fn.c_str(),"mps");
    ClpSimplex solution;
    solution.loadProblem(*m.getMatrixByCol(),m.getColLower(),m.getColUpper(),
			 m.getObjCoefficients(),
			 m.getRowLower(),m.getRowUpper());
    solution.setDblParam(ClpObjOffset,m.objectiveOffset());
    solution.dual();
    // get quadratic part
    int * start=NULL;
    int * column = NULL;
    double * element = NULL;
    m.readQuadraticMps(NULL,start,column,element,2);
    // Load up objective
    solution.loadQuadraticObjective(solution.numberColumns(),start,column,element);
    delete [] start;
    delete [] column;
    delete [] element;
    solution.primal(1);
    solution.nonlinearSLP(50,1.0e-4);
    double objValue = solution.getObjValue();
    CoinRelFltEq eq(1.0e-4);
    assert(eq(objValue,0.5));
    solution.primal();
    objValue = solution.getObjValue();
    assert(eq(objValue,0.5));
  }
#endif  
}