Ignore:
Timestamp:
Jan 6, 2019 2:43:06 PM (4 months ago)
Author:
unxusr
Message:

formatting

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/Clp/src/unitTest.cpp

    r2322 r2385  
    7272//#############################################################################
    7373
    74 
    7574// Function Prototypes. Function definitions is in this file.
    76 void testingMessage( const char * const msg );
     75void testingMessage(const char *const msg);
    7776#if defined(COIN_HAS_AMD) || defined(COIN_HAS_CHOLMOD) || defined(COIN_HAS_GLPK)
    7877static int barrierAvailable = 1;
     
    9089#define NUMBER_ALGORITHMS 12
    9190// If you just want a subset then set some to 1
    92 static int switchOff[NUMBER_ALGORITHMS] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
     91static int switchOff[NUMBER_ALGORITHMS] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 };
    9392// shortName - 0 no , 1 yes
    94 ClpSolve setupForSolve(int algorithm, std::string & nameAlgorithm,
    95                        int shortName)
     93ClpSolve setupForSolve(int algorithm, std::string &nameAlgorithm,
     94  int shortName)
    9695{
    97      ClpSolve solveOptions;
    98      /* algorithms are
     96  ClpSolve solveOptions;
     97  /* algorithms are
    9998        0 barrier
    10099        1 dual with volumne crash
     
    105104        10,11 primal with 70, dual with volume
    106105     */
    107      switch (algorithm) {
    108      case 0:
    109           if (shortName)
    110                nameAlgorithm = "ba";
    111           else
    112                nameAlgorithm = "nameBarrier";
    113           solveOptions.setSolveType(ClpSolve::useBarrier);
    114           if (barrierAvailable == 1)
    115                solveOptions.setSpecialOption(4, 4);
    116           else if (barrierAvailable == 2)
    117                solveOptions.setSpecialOption(4, 2);
    118           break;
    119      case 1:
     106  switch (algorithm) {
     107  case 0:
     108    if (shortName)
     109      nameAlgorithm = "ba";
     110    else
     111      nameAlgorithm = "nameBarrier";
     112    solveOptions.setSolveType(ClpSolve::useBarrier);
     113    if (barrierAvailable == 1)
     114      solveOptions.setSpecialOption(4, 4);
     115    else if (barrierAvailable == 2)
     116      solveOptions.setSpecialOption(4, 2);
     117    break;
     118  case 1:
    120119#ifdef COIN_HAS_VOL
    121           if (shortName)
    122                nameAlgorithm = "du-vol-50";
    123           else
    124                nameAlgorithm = "dual-volume-50";
    125           solveOptions.setSolveType(ClpSolve::useDual);
    126           solveOptions.setSpecialOption(0, 2, 50); // volume
     120    if (shortName)
     121      nameAlgorithm = "du-vol-50";
     122    else
     123      nameAlgorithm = "dual-volume-50";
     124    solveOptions.setSolveType(ClpSolve::useDual);
     125    solveOptions.setSpecialOption(0, 2, 50); // volume
    127126#else
    128           solveOptions.setSolveType(ClpSolve::notImplemented);
    129 #endif
    130           break;
    131      case 2:
    132           if (shortName)
    133                nameAlgorithm = "du-cr";
    134           else
    135                nameAlgorithm = "dual-crash";
    136           solveOptions.setSolveType(ClpSolve::useDual);
    137           solveOptions.setSpecialOption(0, 1);
    138           break;
    139      case 3:
    140           if (shortName)
    141                nameAlgorithm = "du";
    142           else
    143                nameAlgorithm = "dual";
    144           solveOptions.setSolveType(ClpSolve::useDual);
    145           break;
    146      case 4:
    147           if (shortName)
    148                nameAlgorithm = "pr-cr";
    149           else
    150                nameAlgorithm = "primal-crash";
    151           solveOptions.setSolveType(ClpSolve::usePrimal);
    152           solveOptions.setSpecialOption(1, 1);
    153           break;
    154      case 5:
    155           if (shortName)
    156                nameAlgorithm = "pr";
    157           else
    158                nameAlgorithm = "primal";
    159           solveOptions.setSolveType(ClpSolve::usePrimal);
    160           break;
    161      case 6:
    162           if (shortName)
    163                nameAlgorithm = "au-cr";
    164           else
    165                nameAlgorithm = "either-crash";
    166           solveOptions.setSolveType(ClpSolve::automatic);
    167           solveOptions.setSpecialOption(1, 1);
    168           break;
    169      case 7:
    170           if (shortName)
    171                nameAlgorithm = "au";
    172           else
    173                nameAlgorithm = "either";
    174           solveOptions.setSolveType(ClpSolve::automatic);
    175           break;
    176      case 8:
    177           if (shortName)
    178                nameAlgorithm = "pr-id-1";
    179           else
    180                nameAlgorithm = "primal-idiot-1";
    181           solveOptions.setSolveType(ClpSolve::usePrimalorSprint);
    182           solveOptions.setSpecialOption(1, 2, 1); // idiot
    183           break;
    184      case 9:
    185           if (shortName)
    186                nameAlgorithm = "pr-id-5";
    187           else
    188                nameAlgorithm = "primal-idiot-5";
    189           solveOptions.setSolveType(ClpSolve::usePrimalorSprint);
    190           solveOptions.setSpecialOption(1, 2, 5); // idiot
    191           break;
    192      case 10:
    193           if (shortName)
    194                nameAlgorithm = "pr-id-70";
    195           else
    196                nameAlgorithm = "primal-idiot-70";
    197           solveOptions.setSolveType(ClpSolve::usePrimalorSprint);
    198           solveOptions.setSpecialOption(1, 2, 70); // idiot
    199           break;
    200      case 11:
     127    solveOptions.setSolveType(ClpSolve::notImplemented);
     128#endif
     129    break;
     130  case 2:
     131    if (shortName)
     132      nameAlgorithm = "du-cr";
     133    else
     134      nameAlgorithm = "dual-crash";
     135    solveOptions.setSolveType(ClpSolve::useDual);
     136    solveOptions.setSpecialOption(0, 1);
     137    break;
     138  case 3:
     139    if (shortName)
     140      nameAlgorithm = "du";
     141    else
     142      nameAlgorithm = "dual";
     143    solveOptions.setSolveType(ClpSolve::useDual);
     144    break;
     145  case 4:
     146    if (shortName)
     147      nameAlgorithm = "pr-cr";
     148    else
     149      nameAlgorithm = "primal-crash";
     150    solveOptions.setSolveType(ClpSolve::usePrimal);
     151    solveOptions.setSpecialOption(1, 1);
     152    break;
     153  case 5:
     154    if (shortName)
     155      nameAlgorithm = "pr";
     156    else
     157      nameAlgorithm = "primal";
     158    solveOptions.setSolveType(ClpSolve::usePrimal);
     159    break;
     160  case 6:
     161    if (shortName)
     162      nameAlgorithm = "au-cr";
     163    else
     164      nameAlgorithm = "either-crash";
     165    solveOptions.setSolveType(ClpSolve::automatic);
     166    solveOptions.setSpecialOption(1, 1);
     167    break;
     168  case 7:
     169    if (shortName)
     170      nameAlgorithm = "au";
     171    else
     172      nameAlgorithm = "either";
     173    solveOptions.setSolveType(ClpSolve::automatic);
     174    break;
     175  case 8:
     176    if (shortName)
     177      nameAlgorithm = "pr-id-1";
     178    else
     179      nameAlgorithm = "primal-idiot-1";
     180    solveOptions.setSolveType(ClpSolve::usePrimalorSprint);
     181    solveOptions.setSpecialOption(1, 2, 1); // idiot
     182    break;
     183  case 9:
     184    if (shortName)
     185      nameAlgorithm = "pr-id-5";
     186    else
     187      nameAlgorithm = "primal-idiot-5";
     188    solveOptions.setSolveType(ClpSolve::usePrimalorSprint);
     189    solveOptions.setSpecialOption(1, 2, 5); // idiot
     190    break;
     191  case 10:
     192    if (shortName)
     193      nameAlgorithm = "pr-id-70";
     194    else
     195      nameAlgorithm = "primal-idiot-70";
     196    solveOptions.setSolveType(ClpSolve::usePrimalorSprint);
     197    solveOptions.setSpecialOption(1, 2, 70); // idiot
     198    break;
     199  case 11:
    201200#ifdef COIN_HAS_VOL
    202           if (shortName)
    203                nameAlgorithm = "du-vol";
    204           else
    205                nameAlgorithm = "dual-volume";
    206           solveOptions.setSolveType(ClpSolve::useDual);
    207           solveOptions.setSpecialOption(0, 2, 3000); // volume
     201    if (shortName)
     202      nameAlgorithm = "du-vol";
     203    else
     204      nameAlgorithm = "dual-volume";
     205    solveOptions.setSolveType(ClpSolve::useDual);
     206    solveOptions.setSpecialOption(0, 2, 3000); // volume
    208207#else
    209           solveOptions.setSolveType(ClpSolve::notImplemented);
    210 #endif
    211           break;
    212      default:
    213           abort();
    214      }
    215      if (shortName) {
    216           // can switch off
    217           if (switchOff[algorithm])
    218                solveOptions.setSolveType(ClpSolve::notImplemented);
    219      }
    220      return solveOptions;
     208    solveOptions.setSolveType(ClpSolve::notImplemented);
     209#endif
     210    break;
     211  default:
     212    abort();
     213  }
     214  if (shortName) {
     215    // can switch off
     216    if (switchOff[algorithm])
     217      solveOptions.setSolveType(ClpSolve::notImplemented);
     218  }
     219  return solveOptions;
    221220}
    222 static void printSol(ClpSimplex & model)
     221static void printSol(ClpSimplex &model)
    223222{
    224      int numberRows = model.numberRows();
    225      int numberColumns = model.numberColumns();
    226 
    227      double * rowPrimal = model.primalRowSolution();
    228      double * rowDual = model.dualRowSolution();
    229      double * rowLower = model.rowLower();
    230      double * rowUpper = model.rowUpper();
    231 
    232      int iRow;
    233      double objValue = model.getObjValue();
    234      printf("Objvalue %g Rows (%d)\n", objValue, numberRows);
    235      for (iRow = 0; iRow < numberRows; iRow++) {
    236           printf("%d primal %g dual %g low %g up %g\n",
    237                  iRow, rowPrimal[iRow], rowDual[iRow],
    238                  rowLower[iRow], rowUpper[iRow]);
    239      }
    240      double * columnPrimal = model.primalColumnSolution();
    241      double * columnDual = model.dualColumnSolution();
    242      double * columnLower = model.columnLower();
    243      double * columnUpper = model.columnUpper();
    244      double offset;
    245      //const double * gradient = model.objectiveAsObject()->gradient(&model,
    246      //                                                       columnPrimal,offset,true,1);
    247      const double * gradient = model.objective(columnPrimal, offset);
    248      int iColumn;
    249      objValue = -offset - model.objectiveOffset();
    250      printf("offset %g (%g)\n", offset, model.objectiveOffset());
    251      printf("Columns (%d)\n", numberColumns);
    252      for (iColumn = 0; iColumn < numberColumns; iColumn++) {
    253           printf("%d primal %g dual %g low %g up %g\n",
    254                  iColumn, columnPrimal[iColumn], columnDual[iColumn],
    255                  columnLower[iColumn], columnUpper[iColumn]);
    256           objValue += columnPrimal[iColumn] * gradient[iColumn];
    257           if (fabs(columnPrimal[iColumn]*gradient[iColumn]) > 1.0e-8)
    258                printf("obj -> %g gradient %g\n", objValue, gradient[iColumn]);
    259      }
    260      printf("Computed objective %g\n", objValue);
     223  int numberRows = model.numberRows();
     224  int numberColumns = model.numberColumns();
     225
     226  double *rowPrimal = model.primalRowSolution();
     227  double *rowDual = model.dualRowSolution();
     228  double *rowLower = model.rowLower();
     229  double *rowUpper = model.rowUpper();
     230
     231  int iRow;
     232  double objValue = model.getObjValue();
     233  printf("Objvalue %g Rows (%d)\n", objValue, numberRows);
     234  for (iRow = 0; iRow < numberRows; iRow++) {
     235    printf("%d primal %g dual %g low %g up %g\n",
     236      iRow, rowPrimal[iRow], rowDual[iRow],
     237      rowLower[iRow], rowUpper[iRow]);
     238  }
     239  double *columnPrimal = model.primalColumnSolution();
     240  double *columnDual = model.dualColumnSolution();
     241  double *columnLower = model.columnLower();
     242  double *columnUpper = model.columnUpper();
     243  double offset;
     244  //const double * gradient = model.objectiveAsObject()->gradient(&model,
     245  //                                                       columnPrimal,offset,true,1);
     246  const double *gradient = model.objective(columnPrimal, offset);
     247  int iColumn;
     248  objValue = -offset - model.objectiveOffset();
     249  printf("offset %g (%g)\n", offset, model.objectiveOffset());
     250  printf("Columns (%d)\n", numberColumns);
     251  for (iColumn = 0; iColumn < numberColumns; iColumn++) {
     252    printf("%d primal %g dual %g low %g up %g\n",
     253      iColumn, columnPrimal[iColumn], columnDual[iColumn],
     254      columnLower[iColumn], columnUpper[iColumn]);
     255    objValue += columnPrimal[iColumn] * gradient[iColumn];
     256    if (fabs(columnPrimal[iColumn] * gradient[iColumn]) > 1.0e-8)
     257      printf("obj -> %g gradient %g\n", objValue, gradient[iColumn]);
     258  }
     259  printf("Computed objective %g\n", objValue);
    261260}
    262261
    263 void usage(const std::string& key)
     262void usage(const std::string &key)
    264263{
    265      std::cerr
    266                << "Undefined parameter \"" << key << "\".\n"
    267                << "Correct usage: \n"
    268                << "  clp [-dirSample=V1] [-dirNetlib=V2] [-netlib]\n"
    269                << "  where:\n"
    270                << "    -dirSample: directory containing mps test files\n"
    271                << "        Default value V1=\"../../Data/Sample\"\n"
    272                << "    -dirNetlib: directory containing netlib files\"\n"
    273                << "        Default value V2=\"../../Data/Netlib\"\n"
    274                << "    -netlib\n"
    275                << "        If specified, then netlib testset run as well as the nitTest.\n";
     264  std::cerr
     265    << "Undefined parameter \"" << key << "\".\n"
     266    << "Correct usage: \n"
     267    << "  clp [-dirSample=V1] [-dirNetlib=V2] [-netlib]\n"
     268    << "  where:\n"
     269    << "    -dirSample: directory containing mps test files\n"
     270    << "        Default value V1=\"../../Data/Sample\"\n"
     271    << "    -dirNetlib: directory containing netlib files\"\n"
     272    << "        Default value V2=\"../../Data/Netlib\"\n"
     273    << "    -netlib\n"
     274    << "        If specified, then netlib testset run as well as the nitTest.\n";
    276275}
    277276#if FACTORIZATION_STATISTICS
    278 int loSizeX=-1;
    279 int hiSizeX=1000000;
     277int loSizeX = -1;
     278int hiSizeX = 1000000;
    280279#endif
    281280//----------------------------------------------------------------
     
    286285#define AnySimplex AbcSimplex
    287286#endif
    288 int mainTest (int argc, const char *argv[], int algorithm,
    289               AnySimplex empty, ClpSolve solveOptionsIn,
    290               int switchOffValue, bool doVector)
     287int mainTest(int argc, const char *argv[], int algorithm,
     288  AnySimplex empty, ClpSolve solveOptionsIn,
     289  int switchOffValue, bool doVector)
    291290{
    292      int i;
    293 
    294      if (switchOffValue > 0) {
    295           // switch off some
    296           int iTest;
    297           for (iTest = 0; iTest < NUMBER_ALGORITHMS; iTest++) {
     291  int i;
     292
     293  if (switchOffValue > 0) {
     294    // switch off some
     295    int iTest;
     296    for (iTest = 0; iTest < NUMBER_ALGORITHMS; iTest++) {
    298297#ifndef NDEBUG
    299                int bottom = switchOffValue % 10;
    300                assert (bottom == 0 || bottom == 1);
    301 #endif
    302                switchOffValue /= 10;
    303                switchOff[iTest] = 0;
     298      int bottom = switchOffValue % 10;
     299      assert(bottom == 0 || bottom == 1);
     300#endif
     301      switchOffValue /= 10;
     302      switchOff[iTest] = 0;
     303    }
     304  }
     305  int numberFailures = 0;
     306  // define valid parameter keywords
     307  std::set< std::string > definedKeyWords;
     308  definedKeyWords.insert("-dirSample");
     309  definedKeyWords.insert("-dirNetlib");
     310  definedKeyWords.insert("-netlib");
     311
     312  // Create a map of parameter keys and associated data
     313  std::map< std::string, std::string > parms;
     314  for (i = 1; i < argc; i++) {
     315    std::string parm(argv[i]);
     316    std::string key, value;
     317    std::string::size_type eqPos = parm.find('=');
     318
     319    // Does parm contain and '='
     320    if (eqPos == std::string::npos) {
     321      //Parm does not contain '='
     322      key = parm;
     323    } else {
     324      key = parm.substr(0, eqPos);
     325      value = parm.substr(eqPos + 1);
     326    }
     327
     328    // Is specifed key valid?
     329    if (definedKeyWords.find(key) == definedKeyWords.end()) {
     330      // invalid key word.
     331      // Write help text
     332      usage(key);
     333      return 1;
     334    }
     335    parms[key] = value;
     336  }
     337
     338  const char dirsep = CoinFindDirSeparator();
     339  // Set directory containing mps data files.
     340  std::string dirSample;
     341  if (parms.find("-dirSample") != parms.end())
     342    dirSample = parms["-dirSample"];
     343  else
     344    dirSample = dirsep == '/' ? "../../Data/Sample/" : "..\\..\\Data\\Sample\\";
     345
     346  // Set directory containing netlib data files.
     347  std::string dirNetlib;
     348  if (parms.find("-dirNetlib") != parms.end())
     349    dirNetlib = parms["-dirNetlib"];
     350  else
     351    dirNetlib = dirsep == '/' ? "../../Data/Netlib/" : "..\\..\\Data\\Netlib\\";
     352#if FACTORIZATION_STATISTICS == 0
     353  if (!empty.numberRows()) {
     354    testingMessage("Testing ClpSimplex\n");
     355    ClpSimplexUnitTest(dirSample);
     356  }
     357#endif
     358  if (parms.find("-netlib") != parms.end() || empty.numberRows()) {
     359    unsigned int m;
     360    std::string sizeLoHi;
     361#if FACTORIZATION_STATISTICS
     362    double ftranTwiddleFactor1XChoice[3] = { 0.9, 1.0, 1.1 };
     363    double ftranTwiddleFactor2XChoice[3] = { 0.9, 1.0, 1.1 };
     364    double ftranFTTwiddleFactor1XChoice[3] = { 0.9, 1.0, 1.1 };
     365    double ftranFTTwiddleFactor2XChoice[3] = { 0.9, 1.0, 1.1 };
     366    double btranTwiddleFactor1XChoice[3] = { 0.9, 1.0, 1.1 };
     367    double btranTwiddleFactor2XChoice[3] = { 0.9, 1.0, 1.1 };
     368    double ftranFullTwiddleFactor1XChoice[3] = { 0.9, 1.0, 1.1 };
     369    double ftranFullTwiddleFactor2XChoice[3] = { 0.9, 1.0, 1.1 };
     370    double btranFullTwiddleFactor1XChoice[3] = { 0.9, 1.0, 1.1 };
     371    double btranFullTwiddleFactor2XChoice[3] = { 0.9, 1.0, 1.1 };
     372    double denseThresholdXChoice[3] = { 2, 20, 40 };
     373    double twoThresholdXChoice[3] = { 800, 1000, 1200 };
     374    double minRowsSparseChoice[3] = { 250, 300, 350 };
     375    double largeRowsSparseChoice[3] = { 8000, 10000, 12000 };
     376    double mediumRowsDividerChoice[3] = { 5, 6, 7 };
     377    double mediumRowsMinCountChoice[3] = { 300, 500, 600 };
     378    double largeRowsCountChoice[3] = { 700, 1000, 1300 };
     379    double times[3 * 3 * 3 * 3 * 3];
     380    memset(times, 0, sizeof(times));
     381#define whichParam(za, zname)            \
     382  const char *choice##za = #zname;       \
     383  const double *use##za = zname##Choice; \
     384  double *external##za = &zname
     385    whichParam(A, denseThresholdX);
     386    whichParam(B, twoThresholdX);
     387    whichParam(C, minRowsSparse);
     388    whichParam(D, mediumRowsDivider);
     389    whichParam(E, mediumRowsMinCount);
     390#endif
     391    // Define test problems:
     392    //   mps names,
     393    //   maximization or minimization,
     394    //   Number of rows and columns in problem, and
     395    //   objective function value
     396    std::vector< std::string > mpsName;
     397    std::vector< bool > min;
     398    std::vector< int > nRows;
     399    std::vector< int > nCols;
     400    std::vector< double > objValue;
     401    std::vector< double > objValueTol;
     402    // 100 added means no presolve
     403    std::vector< int > bestStrategy;
     404    if (empty.numberRows()) {
     405      std::string alg;
     406      for (int iTest = 0; iTest < NUMBER_ALGORITHMS; iTest++) {
     407        ClpSolve solveOptions = setupForSolve(iTest, alg, 0);
     408        printf("%d %s ", iTest, alg.c_str());
     409        if (switchOff[iTest])
     410          printf("skipped by user\n");
     411        else if (solveOptions.getSolveType() == ClpSolve::notImplemented)
     412          printf("skipped as not available\n");
     413        else
     414          printf("will be tested\n");
     415      }
     416    }
     417    if (!empty.numberRows()) {
     418      mpsName.push_back("25fv47");
     419      min.push_back(true);
     420      nRows.push_back(822);
     421      nCols.push_back(1571);
     422      objValueTol.push_back(1.e-8);
     423      objValue.push_back(5.5018458883E+03);
     424      bestStrategy.push_back(0);
     425      mpsName.push_back("80bau3b");
     426      min.push_back(true);
     427      nRows.push_back(2263);
     428      nCols.push_back(9799);
     429      objValueTol.push_back(1.e-8);
     430      objValue.push_back(9.8722419241E+05);
     431      bestStrategy.push_back(3);
     432      mpsName.push_back("blend");
     433      min.push_back(true);
     434      nRows.push_back(75);
     435      nCols.push_back(83);
     436      objValueTol.push_back(1.e-8);
     437      objValue.push_back(-3.0812149846e+01);
     438      bestStrategy.push_back(3);
     439      mpsName.push_back("pilotnov");
     440      min.push_back(true);
     441      nRows.push_back(976);
     442      nCols.push_back(2172);
     443      objValueTol.push_back(1.e-8);
     444      objValue.push_back(-4.4972761882e+03);
     445      bestStrategy.push_back(3);
     446      mpsName.push_back("maros-r7");
     447      min.push_back(true);
     448      nRows.push_back(3137);
     449      nCols.push_back(9408);
     450      objValueTol.push_back(1.e-8);
     451      objValue.push_back(1.4971851665e+06);
     452      bestStrategy.push_back(2);
     453      mpsName.push_back("pilot");
     454      min.push_back(true);
     455      nRows.push_back(1442);
     456      nCols.push_back(3652);
     457      objValueTol.push_back(1.e-5);
     458      objValue.push_back(/*-5.5740430007e+02*/ -557.48972927292);
     459      bestStrategy.push_back(3);
     460      mpsName.push_back("pilot4");
     461      min.push_back(true);
     462      nRows.push_back(411);
     463      nCols.push_back(1000);
     464      objValueTol.push_back(5.e-5);
     465      objValue.push_back(-2.5811392641e+03);
     466      bestStrategy.push_back(3);
     467      mpsName.push_back("pilot87");
     468      min.push_back(true);
     469      nRows.push_back(2031);
     470      nCols.push_back(4883);
     471      objValueTol.push_back(1.e-4);
     472      objValue.push_back(3.0171072827e+02);
     473      bestStrategy.push_back(0);
     474      mpsName.push_back("adlittle");
     475      min.push_back(true);
     476      nRows.push_back(57);
     477      nCols.push_back(97);
     478      objValueTol.push_back(1.e-8);
     479      objValue.push_back(2.2549496316e+05);
     480      bestStrategy.push_back(3);
     481      mpsName.push_back("afiro");
     482      min.push_back(true);
     483      nRows.push_back(28);
     484      nCols.push_back(32);
     485      objValueTol.push_back(1.e-8);
     486      objValue.push_back(-4.6475314286e+02);
     487      bestStrategy.push_back(3);
     488      mpsName.push_back("agg");
     489      min.push_back(true);
     490      nRows.push_back(489);
     491      nCols.push_back(163);
     492      objValueTol.push_back(1.e-8);
     493      objValue.push_back(-3.5991767287e+07);
     494      bestStrategy.push_back(3);
     495      mpsName.push_back("agg2");
     496      min.push_back(true);
     497      nRows.push_back(517);
     498      nCols.push_back(302);
     499      objValueTol.push_back(1.e-8);
     500      objValue.push_back(-2.0239252356e+07);
     501      bestStrategy.push_back(3);
     502      mpsName.push_back("agg3");
     503      min.push_back(true);
     504      nRows.push_back(517);
     505      nCols.push_back(302);
     506      objValueTol.push_back(1.e-8);
     507      objValue.push_back(1.0312115935e+07);
     508      bestStrategy.push_back(4);
     509      mpsName.push_back("bandm");
     510      min.push_back(true);
     511      nRows.push_back(306);
     512      nCols.push_back(472);
     513      objValueTol.push_back(1.e-8);
     514      objValue.push_back(-1.5862801845e+02);
     515      bestStrategy.push_back(2);
     516      mpsName.push_back("beaconfd");
     517      min.push_back(true);
     518      nRows.push_back(174);
     519      nCols.push_back(262);
     520      objValueTol.push_back(1.e-8);
     521      objValue.push_back(3.3592485807e+04);
     522      bestStrategy.push_back(0);
     523      mpsName.push_back("bnl1");
     524      min.push_back(true);
     525      nRows.push_back(644);
     526      nCols.push_back(1175);
     527      objValueTol.push_back(1.e-8);
     528      objValue.push_back(1.9776295615E+03);
     529      bestStrategy.push_back(3);
     530      mpsName.push_back("bnl2");
     531      min.push_back(true);
     532      nRows.push_back(2325);
     533      nCols.push_back(3489);
     534      objValueTol.push_back(1.e-8);
     535      objValue.push_back(1.8112365404e+03);
     536      bestStrategy.push_back(3);
     537      mpsName.push_back("boeing1");
     538      min.push_back(true);
     539      nRows.push_back(/*351*/ 352);
     540      nCols.push_back(384);
     541      objValueTol.push_back(1.e-8);
     542      objValue.push_back(-3.3521356751e+02);
     543      bestStrategy.push_back(3);
     544      mpsName.push_back("boeing2");
     545      min.push_back(true);
     546      nRows.push_back(167);
     547      nCols.push_back(143);
     548      objValueTol.push_back(1.e-8);
     549      objValue.push_back(-3.1501872802e+02);
     550      bestStrategy.push_back(3);
     551      mpsName.push_back("bore3d");
     552      min.push_back(true);
     553      nRows.push_back(234);
     554      nCols.push_back(315);
     555      objValueTol.push_back(1.e-8);
     556      objValue.push_back(1.3730803942e+03);
     557      bestStrategy.push_back(3);
     558      mpsName.push_back("brandy");
     559      min.push_back(true);
     560      nRows.push_back(221);
     561      nCols.push_back(249);
     562      objValueTol.push_back(1.e-8);
     563      objValue.push_back(1.5185098965e+03);
     564      bestStrategy.push_back(3);
     565      mpsName.push_back("capri");
     566      min.push_back(true);
     567      nRows.push_back(272);
     568      nCols.push_back(353);
     569      objValueTol.push_back(1.e-8);
     570      objValue.push_back(2.6900129138e+03);
     571      bestStrategy.push_back(3);
     572      mpsName.push_back("cycle");
     573      min.push_back(true);
     574      nRows.push_back(1904);
     575      nCols.push_back(2857);
     576      objValueTol.push_back(1.e-9);
     577      objValue.push_back(-5.2263930249e+00);
     578      bestStrategy.push_back(3);
     579      mpsName.push_back("czprob");
     580      min.push_back(true);
     581      nRows.push_back(930);
     582      nCols.push_back(3523);
     583      objValueTol.push_back(1.e-8);
     584      objValue.push_back(2.1851966989e+06);
     585      bestStrategy.push_back(3);
     586      mpsName.push_back("d2q06c");
     587      min.push_back(true);
     588      nRows.push_back(2172);
     589      nCols.push_back(5167);
     590      objValueTol.push_back(1.e-7);
     591      objValue.push_back(122784.21557456);
     592      bestStrategy.push_back(0);
     593      mpsName.push_back("d6cube");
     594      min.push_back(true);
     595      nRows.push_back(416);
     596      nCols.push_back(6184);
     597      objValueTol.push_back(1.e-7);
     598      objValue.push_back(3.1549166667e+02);
     599      bestStrategy.push_back(3);
     600      mpsName.push_back("degen2");
     601      min.push_back(true);
     602      nRows.push_back(445);
     603      nCols.push_back(534);
     604      objValueTol.push_back(1.e-8);
     605      objValue.push_back(-1.4351780000e+03);
     606      bestStrategy.push_back(3);
     607      mpsName.push_back("degen3");
     608      min.push_back(true);
     609      nRows.push_back(1504);
     610      nCols.push_back(1818);
     611      objValueTol.push_back(1.e-8);
     612      objValue.push_back(-9.8729400000e+02);
     613      bestStrategy.push_back(2);
     614      mpsName.push_back("dfl001");
     615      min.push_back(true);
     616      nRows.push_back(6072);
     617      nCols.push_back(12230);
     618      objValueTol.push_back(1.e-5);
     619      objValue.push_back(1.1266396047E+07);
     620      bestStrategy.push_back(5);
     621      mpsName.push_back("e226");
     622      min.push_back(true);
     623      nRows.push_back(224);
     624      nCols.push_back(282);
     625      objValueTol.push_back(1.e-8);
     626      objValue.push_back(-1.8751929066e+01 + 7.113);
     627      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.
     628      mpsName.push_back("etamacro");
     629      min.push_back(true);
     630      nRows.push_back(401);
     631      nCols.push_back(688);
     632      objValueTol.push_back(1.e-6);
     633      objValue.push_back(-7.5571521774e+02);
     634      bestStrategy.push_back(3);
     635      mpsName.push_back("fffff800");
     636      min.push_back(true);
     637      nRows.push_back(525);
     638      nCols.push_back(854);
     639      objValueTol.push_back(1.e-6);
     640      objValue.push_back(5.5567961165e+05);
     641      bestStrategy.push_back(3);
     642      mpsName.push_back("finnis");
     643      min.push_back(true);
     644      nRows.push_back(498);
     645      nCols.push_back(614);
     646      objValueTol.push_back(1.e-6);
     647      objValue.push_back(1.7279096547e+05);
     648      bestStrategy.push_back(3);
     649      mpsName.push_back("fit1d");
     650      min.push_back(true);
     651      nRows.push_back(25);
     652      nCols.push_back(1026);
     653      objValueTol.push_back(1.e-8);
     654      objValue.push_back(-9.1463780924e+03);
     655      bestStrategy.push_back(3 + 100);
     656      mpsName.push_back("fit1p");
     657      min.push_back(true);
     658      nRows.push_back(628);
     659      nCols.push_back(1677);
     660      objValueTol.push_back(1.e-8);
     661      objValue.push_back(9.1463780924e+03);
     662      bestStrategy.push_back(5 + 100);
     663      mpsName.push_back("fit2d");
     664      min.push_back(true);
     665      nRows.push_back(26);
     666      nCols.push_back(10500);
     667      objValueTol.push_back(1.e-8);
     668      objValue.push_back(-6.8464293294e+04);
     669      bestStrategy.push_back(3 + 100);
     670      mpsName.push_back("fit2p");
     671      min.push_back(true);
     672      nRows.push_back(3001);
     673      nCols.push_back(13525);
     674      objValueTol.push_back(1.e-9);
     675      objValue.push_back(6.8464293232e+04);
     676      bestStrategy.push_back(5 + 100);
     677      mpsName.push_back("forplan");
     678      min.push_back(true);
     679      nRows.push_back(162);
     680      nCols.push_back(421);
     681      objValueTol.push_back(1.e-6);
     682      objValue.push_back(-6.6421873953e+02);
     683      bestStrategy.push_back(3);
     684      mpsName.push_back("ganges");
     685      min.push_back(true);
     686      nRows.push_back(1310);
     687      nCols.push_back(1681);
     688      objValueTol.push_back(1.e-5);
     689      objValue.push_back(-1.0958636356e+05);
     690      bestStrategy.push_back(3);
     691      mpsName.push_back("gfrd-pnc");
     692      min.push_back(true);
     693      nRows.push_back(617);
     694      nCols.push_back(1092);
     695      objValueTol.push_back(1.e-8);
     696      objValue.push_back(6.9022359995e+06);
     697      bestStrategy.push_back(3);
     698      mpsName.push_back("greenbea");
     699      min.push_back(true);
     700      nRows.push_back(2393);
     701      nCols.push_back(5405);
     702      objValueTol.push_back(1.e-8);
     703      objValue.push_back(/*-7.2462405908e+07*/ -72555248.129846);
     704      bestStrategy.push_back(3);
     705      mpsName.push_back("greenbeb");
     706      min.push_back(true);
     707      nRows.push_back(2393);
     708      nCols.push_back(5405);
     709      objValueTol.push_back(1.e-8);
     710      objValue.push_back(/*-4.3021476065e+06*/ -4302260.2612066);
     711      bestStrategy.push_back(3);
     712      mpsName.push_back("grow15");
     713      min.push_back(true);
     714      nRows.push_back(301);
     715      nCols.push_back(645);
     716      objValueTol.push_back(1.e-8);
     717      objValue.push_back(-1.0687094129e+08);
     718      bestStrategy.push_back(4 + 100);
     719      mpsName.push_back("grow22");
     720      min.push_back(true);
     721      nRows.push_back(441);
     722      nCols.push_back(946);
     723      objValueTol.push_back(1.e-8);
     724      objValue.push_back(-1.6083433648e+08);
     725      bestStrategy.push_back(4 + 100);
     726      mpsName.push_back("grow7");
     727      min.push_back(true);
     728      nRows.push_back(141);
     729      nCols.push_back(301);
     730      objValueTol.push_back(1.e-8);
     731      objValue.push_back(-4.7787811815e+07);
     732      bestStrategy.push_back(4 + 100);
     733      mpsName.push_back("israel");
     734      min.push_back(true);
     735      nRows.push_back(175);
     736      nCols.push_back(142);
     737      objValueTol.push_back(1.e-8);
     738      objValue.push_back(-8.9664482186e+05);
     739      bestStrategy.push_back(2);
     740      mpsName.push_back("kb2");
     741      min.push_back(true);
     742      nRows.push_back(44);
     743      nCols.push_back(41);
     744      objValueTol.push_back(1.e-8);
     745      objValue.push_back(-1.7499001299e+03);
     746      bestStrategy.push_back(3);
     747      mpsName.push_back("lotfi");
     748      min.push_back(true);
     749      nRows.push_back(154);
     750      nCols.push_back(308);
     751      objValueTol.push_back(1.e-8);
     752      objValue.push_back(-2.5264706062e+01);
     753      bestStrategy.push_back(3);
     754      mpsName.push_back("maros");
     755      min.push_back(true);
     756      nRows.push_back(847);
     757      nCols.push_back(1443);
     758      objValueTol.push_back(1.e-8);
     759      objValue.push_back(-5.8063743701e+04);
     760      bestStrategy.push_back(3);
     761      mpsName.push_back("modszk1");
     762      min.push_back(true);
     763      nRows.push_back(688);
     764      nCols.push_back(1620);
     765      objValueTol.push_back(1.e-8);
     766      objValue.push_back(3.2061972906e+02);
     767      bestStrategy.push_back(3);
     768      mpsName.push_back("nesm");
     769      min.push_back(true);
     770      nRows.push_back(663);
     771      nCols.push_back(2923);
     772      objValueTol.push_back(1.e-5);
     773      objValue.push_back(1.4076073035e+07);
     774      bestStrategy.push_back(2);
     775      mpsName.push_back("perold");
     776      min.push_back(true);
     777      nRows.push_back(626);
     778      nCols.push_back(1376);
     779      objValueTol.push_back(1.e-6);
     780      objValue.push_back(-9.3807580773e+03);
     781      bestStrategy.push_back(3);
     782      //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);
     783      //mpsName.push_back("qap15");min.push_back(true);nRows.push_back(6331);nCols.push_back(22275);objValueTol.push_back(1.e-8);objValue.push_back(1.0409940410e+03);bestStrategy.push_back(3);
     784      mpsName.push_back("recipe");
     785      min.push_back(true);
     786      nRows.push_back(92);
     787      nCols.push_back(180);
     788      objValueTol.push_back(1.e-8);
     789      objValue.push_back(-2.6661600000e+02);
     790      bestStrategy.push_back(3);
     791      mpsName.push_back("sc105");
     792      min.push_back(true);
     793      nRows.push_back(106);
     794      nCols.push_back(103);
     795      objValueTol.push_back(1.e-8);
     796      objValue.push_back(-5.2202061212e+01);
     797      bestStrategy.push_back(3);
     798      mpsName.push_back("sc205");
     799      min.push_back(true);
     800      nRows.push_back(206);
     801      nCols.push_back(203);
     802      objValueTol.push_back(1.e-8);
     803      objValue.push_back(-5.2202061212e+01);
     804      bestStrategy.push_back(3);
     805      mpsName.push_back("sc50a");
     806      min.push_back(true);
     807      nRows.push_back(51);
     808      nCols.push_back(48);
     809      objValueTol.push_back(1.e-8);
     810      objValue.push_back(-6.4575077059e+01);
     811      bestStrategy.push_back(3);
     812      mpsName.push_back("sc50b");
     813      min.push_back(true);
     814      nRows.push_back(51);
     815      nCols.push_back(48);
     816      objValueTol.push_back(1.e-8);
     817      objValue.push_back(-7.0000000000e+01);
     818      bestStrategy.push_back(3);
     819      mpsName.push_back("scagr25");
     820      min.push_back(true);
     821      nRows.push_back(472);
     822      nCols.push_back(500);
     823      objValueTol.push_back(1.e-8);
     824      objValue.push_back(-1.4753433061e+07);
     825      bestStrategy.push_back(3);
     826      mpsName.push_back("scagr7");
     827      min.push_back(true);
     828      nRows.push_back(130);
     829      nCols.push_back(140);
     830      objValueTol.push_back(1.e-6);
     831      objValue.push_back(-2.3313892548e+06);
     832      bestStrategy.push_back(3);
     833      mpsName.push_back("scfxm1");
     834      min.push_back(true);
     835      nRows.push_back(331);
     836      nCols.push_back(457);
     837      objValueTol.push_back(1.e-8);
     838      objValue.push_back(1.8416759028e+04);
     839      bestStrategy.push_back(3);
     840      mpsName.push_back("scfxm2");
     841      min.push_back(true);
     842      nRows.push_back(661);
     843      nCols.push_back(914);
     844      objValueTol.push_back(1.e-8);
     845      objValue.push_back(3.6660261565e+04);
     846      bestStrategy.push_back(3);
     847      mpsName.push_back("scfxm3");
     848      min.push_back(true);
     849      nRows.push_back(991);
     850      nCols.push_back(1371);
     851      objValueTol.push_back(1.e-8);
     852      objValue.push_back(5.4901254550e+04);
     853      bestStrategy.push_back(3);
     854      mpsName.push_back("scorpion");
     855      min.push_back(true);
     856      nRows.push_back(389);
     857      nCols.push_back(358);
     858      objValueTol.push_back(1.e-8);
     859      objValue.push_back(1.8781248227e+03);
     860      bestStrategy.push_back(3);
     861      mpsName.push_back("scrs8");
     862      min.push_back(true);
     863      nRows.push_back(491);
     864      nCols.push_back(1169);
     865      objValueTol.push_back(1.e-5);
     866      objValue.push_back(9.0429998619e+02);
     867      bestStrategy.push_back(2);
     868      mpsName.push_back("scsd1");
     869      min.push_back(true);
     870      nRows.push_back(78);
     871      nCols.push_back(760);
     872      objValueTol.push_back(1.e-8);
     873      objValue.push_back(8.6666666743e+00);
     874      bestStrategy.push_back(3 + 100);
     875      mpsName.push_back("scsd6");
     876      min.push_back(true);
     877      nRows.push_back(148);
     878      nCols.push_back(1350);
     879      objValueTol.push_back(1.e-8);
     880      objValue.push_back(5.0500000078e+01);
     881      bestStrategy.push_back(3 + 100);
     882      mpsName.push_back("scsd8");
     883      min.push_back(true);
     884      nRows.push_back(398);
     885      nCols.push_back(2750);
     886      objValueTol.push_back(1.e-7);
     887      objValue.push_back(9.0499999993e+02);
     888      bestStrategy.push_back(1 + 100);
     889      mpsName.push_back("sctap1");
     890      min.push_back(true);
     891      nRows.push_back(301);
     892      nCols.push_back(480);
     893      objValueTol.push_back(1.e-8);
     894      objValue.push_back(1.4122500000e+03);
     895      bestStrategy.push_back(3);
     896      mpsName.push_back("sctap2");
     897      min.push_back(true);
     898      nRows.push_back(1091);
     899      nCols.push_back(1880);
     900      objValueTol.push_back(1.e-8);
     901      objValue.push_back(1.7248071429e+03);
     902      bestStrategy.push_back(3);
     903      mpsName.push_back("sctap3");
     904      min.push_back(true);
     905      nRows.push_back(1481);
     906      nCols.push_back(2480);
     907      objValueTol.push_back(1.e-8);
     908      objValue.push_back(1.4240000000e+03);
     909      bestStrategy.push_back(3);
     910      mpsName.push_back("seba");
     911      min.push_back(true);
     912      nRows.push_back(516);
     913      nCols.push_back(1028);
     914      objValueTol.push_back(1.e-8);
     915      objValue.push_back(1.5711600000e+04);
     916      bestStrategy.push_back(3);
     917      mpsName.push_back("share1b");
     918      min.push_back(true);
     919      nRows.push_back(118);
     920      nCols.push_back(225);
     921      objValueTol.push_back(1.e-8);
     922      objValue.push_back(-7.6589318579e+04);
     923      bestStrategy.push_back(3);
     924      mpsName.push_back("share2b");
     925      min.push_back(true);
     926      nRows.push_back(97);
     927      nCols.push_back(79);
     928      objValueTol.push_back(1.e-8);
     929      objValue.push_back(-4.1573224074e+02);
     930      bestStrategy.push_back(3);
     931      mpsName.push_back("shell");
     932      min.push_back(true);
     933      nRows.push_back(537);
     934      nCols.push_back(1775);
     935      objValueTol.push_back(1.e-8);
     936      objValue.push_back(1.2088253460e+09);
     937      bestStrategy.push_back(3);
     938      mpsName.push_back("ship04l");
     939      min.push_back(true);
     940      nRows.push_back(403);
     941      nCols.push_back(2118);
     942      objValueTol.push_back(1.e-8);
     943      objValue.push_back(1.7933245380e+06);
     944      bestStrategy.push_back(3);
     945      mpsName.push_back("ship04s");
     946      min.push_back(true);
     947      nRows.push_back(403);
     948      nCols.push_back(1458);
     949      objValueTol.push_back(1.e-8);
     950      objValue.push_back(1.7987147004e+06);
     951      bestStrategy.push_back(3);
     952      mpsName.push_back("ship08l");
     953      min.push_back(true);
     954      nRows.push_back(779);
     955      nCols.push_back(4283);
     956      objValueTol.push_back(1.e-8);
     957      objValue.push_back(1.9090552114e+06);
     958      bestStrategy.push_back(3);
     959      mpsName.push_back("ship08s");
     960      min.push_back(true);
     961      nRows.push_back(779);
     962      nCols.push_back(2387);
     963      objValueTol.push_back(1.e-8);
     964      objValue.push_back(1.9200982105e+06);
     965      bestStrategy.push_back(2);
     966      mpsName.push_back("ship12l");
     967      min.push_back(true);
     968      nRows.push_back(1152);
     969      nCols.push_back(5427);
     970      objValueTol.push_back(1.e-8);
     971      objValue.push_back(1.4701879193e+06);
     972      bestStrategy.push_back(3);
     973      mpsName.push_back("ship12s");
     974      min.push_back(true);
     975      nRows.push_back(1152);
     976      nCols.push_back(2763);
     977      objValueTol.push_back(1.e-8);
     978      objValue.push_back(1.4892361344e+06);
     979      bestStrategy.push_back(2);
     980      mpsName.push_back("sierra");
     981      min.push_back(true);
     982      nRows.push_back(1228);
     983      nCols.push_back(2036);
     984      objValueTol.push_back(1.e-8);
     985      objValue.push_back(1.5394362184e+07);
     986      bestStrategy.push_back(3);
     987      mpsName.push_back("stair");
     988      min.push_back(true);
     989      nRows.push_back(357);
     990      nCols.push_back(467);
     991      objValueTol.push_back(1.e-8);
     992      objValue.push_back(-2.5126695119e+02);
     993      bestStrategy.push_back(3);
     994      mpsName.push_back("standata");
     995      min.push_back(true);
     996      nRows.push_back(360);
     997      nCols.push_back(1075);
     998      objValueTol.push_back(1.e-8);
     999      objValue.push_back(1.2576995000e+03);
     1000      bestStrategy.push_back(3);
     1001      //mpsName.push_back("standgub");min.push_back(true);nRows.push_back(362);nCols.push_back(1184);objValueTol.push_back(1.e-8);objValue.push_back(1257.6995); bestStrategy.push_back(3);
     1002      mpsName.push_back("standmps");
     1003      min.push_back(true);
     1004      nRows.push_back(468);
     1005      nCols.push_back(1075);
     1006      objValueTol.push_back(1.e-8);
     1007      objValue.push_back(1.4060175000E+03);
     1008      bestStrategy.push_back(3);
     1009      mpsName.push_back("stocfor1");
     1010      min.push_back(true);
     1011      nRows.push_back(118);
     1012      nCols.push_back(111);
     1013      objValueTol.push_back(1.e-8);
     1014      objValue.push_back(-4.1131976219E+04);
     1015      bestStrategy.push_back(3);
     1016      mpsName.push_back("stocfor2");
     1017      min.push_back(true);
     1018      nRows.push_back(2158);
     1019      nCols.push_back(2031);
     1020      objValueTol.push_back(1.e-8);
     1021      objValue.push_back(-3.9024408538e+04);
     1022      bestStrategy.push_back(3);
     1023      //mpsName.push_back("stocfor3");min.push_back(true);nRows.push_back(16676);nCols.push_back(15695);objValueTol.push_back(1.e-8);objValue.push_back(-3.9976661576e+04);bestStrategy.push_back(3);
     1024      //mpsName.push_back("truss");min.push_back(true);nRows.push_back(1001);nCols.push_back(8806);objValueTol.push_back(1.e-8);objValue.push_back(4.5881584719e+05);bestStrategy.push_back(3);
     1025      mpsName.push_back("tuff");
     1026      min.push_back(true);
     1027      nRows.push_back(334);
     1028      nCols.push_back(587);
     1029      objValueTol.push_back(1.e-8);
     1030      objValue.push_back(2.9214776509e-01);
     1031      bestStrategy.push_back(3);
     1032      mpsName.push_back("vtpbase");
     1033      min.push_back(true);
     1034      nRows.push_back(199);
     1035      nCols.push_back(203);
     1036      objValueTol.push_back(1.e-8);
     1037      objValue.push_back(1.2983146246e+05);
     1038      bestStrategy.push_back(3);
     1039      mpsName.push_back("wood1p");
     1040      min.push_back(true);
     1041      nRows.push_back(245);
     1042      nCols.push_back(2594);
     1043      objValueTol.push_back(5.e-5);
     1044      objValue.push_back(1.4429024116e+00);
     1045      bestStrategy.push_back(3);
     1046      mpsName.push_back("woodw");
     1047      min.push_back(true);
     1048      nRows.push_back(1099);
     1049      nCols.push_back(8405);
     1050      objValueTol.push_back(1.e-8);
     1051      objValue.push_back(1.3044763331E+00);
     1052      bestStrategy.push_back(3);
     1053    } else {
     1054      // Just testing one
     1055      mpsName.push_back(empty.problemName());
     1056      min.push_back(true);
     1057      nRows.push_back(-1);
     1058      nCols.push_back(-1);
     1059      objValueTol.push_back(1.e-8);
     1060      objValue.push_back(0.0);
     1061      bestStrategy.push_back(0);
     1062      int iTest;
     1063      std::string alg;
     1064      for (iTest = 0; iTest < NUMBER_ALGORITHMS; iTest++) {
     1065        ClpSolve solveOptions = setupForSolve(iTest, alg, 0);
     1066        printf("%d %s ", iTest, alg.c_str());
     1067        if (switchOff[iTest])
     1068          printf("skipped by user\n");
     1069        else if (solveOptions.getSolveType() == ClpSolve::notImplemented)
     1070          printf("skipped as not available\n");
     1071        else
     1072          printf("will be tested\n");
     1073      }
     1074    }
     1075
     1076    double timeTaken = 0.0;
     1077    if (!barrierAvailable)
     1078      switchOff[0] = 1;
     1079    // Loop once for each Mps File
     1080    for (m = 0; m < mpsName.size(); m++) {
     1081#if FACTORIZATION_STATISTICS
     1082      if (nRows[m] < loSizeX || nRows[m] >= hiSizeX) {
     1083        std::cerr << "  skipping mps file: " << mpsName[m]
     1084                  << " as " << nRows[m]
     1085                  << " (" << m + 1 << " out of " << mpsName.size() << ")" << std::endl;
     1086        continue;
     1087      }
     1088#endif
     1089      std::cerr << "  processing mps file: " << mpsName[m]
     1090                << " (" << m + 1 << " out of " << mpsName.size() << ")" << std::endl;
     1091      AnySimplex solutionBase = empty;
     1092      std::string fn = dirNetlib + mpsName[m];
     1093      if (!empty.numberRows() || algorithm < 6) {
     1094        // Read data mps file,
     1095        CoinMpsIO mps;
     1096        int nerrors = mps.readMps(fn.c_str(), "mps");
     1097        if (nerrors) {
     1098          std::cerr << "Error " << nerrors << " when reading model from "
     1099                    << fn.c_str() << "! Aborting tests.\n";
     1100          return 1;
     1101        }
     1102        solutionBase.loadProblem(*mps.getMatrixByCol(), mps.getColLower(),
     1103          mps.getColUpper(),
     1104          mps.getObjCoefficients(),
     1105          mps.getRowLower(), mps.getRowUpper());
     1106
     1107        solutionBase.setDblParam(ClpObjOffset, mps.objectiveOffset());
     1108      }
     1109
     1110      // Runs through strategies
     1111      if (algorithm == 6 || algorithm == 7) {
     1112        // algorithms tested are at top of file
     1113        double testTime[NUMBER_ALGORITHMS];
     1114        std::string alg[NUMBER_ALGORITHMS];
     1115        int iTest;
     1116        for (iTest = 0; iTest < NUMBER_ALGORITHMS; iTest++) {
     1117          ClpSolve solveOptions = setupForSolve(iTest, alg[iTest], 1);
     1118          if (solveOptions.getSolveType() != ClpSolve::notImplemented) {
     1119            double time1 = CoinCpuTime();
     1120            ClpSimplex solution = solutionBase;
     1121            if (solution.maximumSeconds() < 0.0)
     1122              solution.setMaximumSeconds(120.0);
     1123            if (doVector) {
     1124              ClpMatrixBase *matrix = solution.clpMatrix();
     1125              if (dynamic_cast< ClpPackedMatrix * >(matrix)) {
     1126                ClpPackedMatrix *clpMatrix = dynamic_cast< ClpPackedMatrix * >(matrix);
     1127                clpMatrix->makeSpecialColumnCopy();
     1128              }
     1129            }
     1130            solution.initialSolve(solveOptions);
     1131            double time2 = CoinCpuTime() - time1;
     1132            testTime[iTest] = time2;
     1133            printf("Finished %s Took %g seconds (%d iterations) - status %d\n",
     1134              mpsName[m].c_str(), time2, solution.problemStatus(), solution.numberIterations());
     1135            if (solution.problemStatus())
     1136              testTime[iTest] = 1.0e20;
     1137          } else {
     1138            testTime[iTest] = 1.0e30;
    3041139          }
    305      }
    306      int numberFailures=0;
    307      // define valid parameter keywords
    308      std::set<std::string> definedKeyWords;
    309      definedKeyWords.insert("-dirSample");
    310      definedKeyWords.insert("-dirNetlib");
    311      definedKeyWords.insert("-netlib");
    312 
    313      // Create a map of parameter keys and associated data
    314      std::map<std::string, std::string> parms;
    315      for ( i = 1; i < argc; i++ ) {
    316           std::string parm(argv[i]);
    317           std::string key, value;
    318           std::string::size_type  eqPos = parm.find('=');
    319 
    320           // Does parm contain and '='
    321           if ( eqPos == std::string::npos ) {
    322                //Parm does not contain '='
    323                key = parm;
    324           } else {
    325                key = parm.substr(0, eqPos);
    326                value = parm.substr(eqPos + 1);
     1140        }
     1141        int iBest = -1;
     1142        double dBest = 1.0e10;
     1143        printf("%s", fn.c_str());
     1144        for (iTest = 0; iTest < NUMBER_ALGORITHMS; iTest++) {
     1145          if (testTime[iTest] < 1.0e30) {
     1146            printf(" %s %g",
     1147              alg[iTest].c_str(), testTime[iTest]);
     1148            if (testTime[iTest] < dBest) {
     1149              dBest = testTime[iTest];
     1150              iBest = iTest;
     1151            }
    3271152          }
    328 
    329           // Is specifed key valid?
    330           if ( definedKeyWords.find(key) == definedKeyWords.end() ) {
    331                // invalid key word.
    332                // Write help text
    333                usage(key);
    334                return 1;
    335           }
    336           parms[key] = value;
    337      }
    338 
    339      const char dirsep =  CoinFindDirSeparator();
    340      // Set directory containing mps data files.
    341      std::string dirSample;
    342      if (parms.find("-dirSample") != parms.end())
    343           dirSample = parms["-dirSample"];
    344      else
    345           dirSample = dirsep == '/' ? "../../Data/Sample/" : "..\\..\\Data\\Sample\\";
    346 
    347      // Set directory containing netlib data files.
    348      std::string dirNetlib;
    349      if (parms.find("-dirNetlib") != parms.end())
    350           dirNetlib = parms["-dirNetlib"];
    351      else
    352           dirNetlib = dirsep == '/' ? "../../Data/Netlib/" : "..\\..\\Data\\Netlib\\";
    353 #if FACTORIZATION_STATISTICS==0
    354      if (!empty.numberRows()) {
    355           testingMessage( "Testing ClpSimplex\n" );
    356           ClpSimplexUnitTest(dirSample);
    357      }
    358 #endif
    359      if (parms.find("-netlib") != parms.end() || empty.numberRows()) {
    360           unsigned int m;
    361           std::string sizeLoHi;
    362 #if FACTORIZATION_STATISTICS
    363           double ftranTwiddleFactor1XChoice[3]={0.9,1.0,1.1};
    364           double ftranTwiddleFactor2XChoice[3]={0.9,1.0,1.1};
    365           double ftranFTTwiddleFactor1XChoice[3]={0.9,1.0,1.1};
    366           double ftranFTTwiddleFactor2XChoice[3]={0.9,1.0,1.1};
    367           double btranTwiddleFactor1XChoice[3]={0.9,1.0,1.1};
    368           double btranTwiddleFactor2XChoice[3]={0.9,1.0,1.1};
    369           double ftranFullTwiddleFactor1XChoice[3]={0.9,1.0,1.1};
    370           double ftranFullTwiddleFactor2XChoice[3]={0.9,1.0,1.1};
    371           double btranFullTwiddleFactor1XChoice[3]={0.9,1.0,1.1};
    372           double btranFullTwiddleFactor2XChoice[3]={0.9,1.0,1.1};
    373           double denseThresholdXChoice[3]={2,20,40};
    374           double twoThresholdXChoice[3]={800,1000,1200};
    375           double minRowsSparseChoice[3]={250,300,350};
    376           double largeRowsSparseChoice[3]={8000,10000,12000};
    377           double mediumRowsDividerChoice[3]={5,6,7};
    378           double mediumRowsMinCountChoice[3]={300,500,600};
    379           double largeRowsCountChoice[3]={700,1000,1300};
    380           double times[3*3*3*3*3];
    381           memset(times,0,sizeof(times));
    382 #define whichParam(za,zname) const char * choice##za=#zname; \
    383        const double * use##za=zname##Choice;                 \
    384        double * external##za=&zname
    385           whichParam(A,denseThresholdX);
    386           whichParam(B,twoThresholdX);
    387           whichParam(C,minRowsSparse);
    388           whichParam(D,mediumRowsDivider);
    389           whichParam(E,mediumRowsMinCount);
    390 #endif
    391           // Define test problems:
    392           //   mps names,
    393           //   maximization or minimization,
    394           //   Number of rows and columns in problem, and
    395           //   objective function value
    396           std::vector<std::string> mpsName;
    397           std::vector<bool> min;
    398           std::vector<int> nRows;
    399           std::vector<int> nCols;
    400           std::vector<double> objValue;
    401           std::vector<double> objValueTol;
    402           // 100 added means no presolve
    403           std::vector<int> bestStrategy;
    404           if(empty.numberRows()) {
    405                std::string alg;
    406                for (int iTest = 0; iTest < NUMBER_ALGORITHMS; iTest++) {
    407                     ClpSolve solveOptions = setupForSolve(iTest, alg, 0);
    408                     printf("%d %s ", iTest, alg.c_str());
    409                     if (switchOff[iTest])
    410                          printf("skipped by user\n");
    411                     else if(solveOptions.getSolveType() == ClpSolve::notImplemented)
    412                          printf("skipped as not available\n");
    413                     else
    414                          printf("will be tested\n");
    415                }
    416           }
    417           if (!empty.numberRows()) {
    418                mpsName.push_back("25fv47");
    419                min.push_back(true);
    420                nRows.push_back(822);
    421                nCols.push_back(1571);
    422                objValueTol.push_back(1.e-8);
    423                objValue.push_back(5.5018458883E+03);
    424                bestStrategy.push_back(0);
    425                mpsName.push_back("80bau3b");
    426                min.push_back(true);
    427                nRows.push_back(2263);
    428                nCols.push_back(9799);
    429                objValueTol.push_back(1.e-8);
    430                objValue.push_back(9.8722419241E+05);
    431                bestStrategy.push_back(3);
    432                mpsName.push_back("blend");
    433                min.push_back(true);
    434                nRows.push_back(75);
    435                nCols.push_back(83);
    436                objValueTol.push_back(1.e-8);
    437                objValue.push_back(-3.0812149846e+01);
    438                bestStrategy.push_back(3);
    439                mpsName.push_back("pilotnov");
    440                min.push_back(true);
    441                nRows.push_back(976);
    442                nCols.push_back(2172);
    443                objValueTol.push_back(1.e-8);
    444                objValue.push_back(-4.4972761882e+03);
    445                bestStrategy.push_back(3);
    446                mpsName.push_back("maros-r7");
    447                min.push_back(true);
    448                nRows.push_back(3137);
    449                nCols.push_back(9408);
    450                objValueTol.push_back(1.e-8);
    451                objValue.push_back(1.4971851665e+06);
    452                bestStrategy.push_back(2);
    453                mpsName.push_back("pilot");
    454                min.push_back(true);
    455                nRows.push_back(1442);
    456                nCols.push_back(3652);
    457                objValueTol.push_back(1.e-5);
    458                objValue.push_back(/*-5.5740430007e+02*/ -557.48972927292);
    459                bestStrategy.push_back(3);
    460                mpsName.push_back("pilot4");
    461                min.push_back(true);
    462                nRows.push_back(411);
    463                nCols.push_back(1000);
    464                objValueTol.push_back(5.e-5);
    465                objValue.push_back(-2.5811392641e+03);
    466                bestStrategy.push_back(3);
    467                mpsName.push_back("pilot87");
    468                min.push_back(true);
    469                nRows.push_back(2031);
    470                nCols.push_back(4883);
    471                objValueTol.push_back(1.e-4);
    472                objValue.push_back(3.0171072827e+02);
    473                bestStrategy.push_back(0);
    474                mpsName.push_back("adlittle");
    475                min.push_back(true);
    476                nRows.push_back(57);
    477                nCols.push_back(97);
    478                objValueTol.push_back(1.e-8);
    479                objValue.push_back(2.2549496316e+05);
    480                bestStrategy.push_back(3);
    481                mpsName.push_back("afiro");
    482                min.push_back(true);
    483                nRows.push_back(28);
    484                nCols.push_back(32);
    485                objValueTol.push_back(1.e-8);
    486                objValue.push_back(-4.6475314286e+02);
    487                bestStrategy.push_back(3);
    488                mpsName.push_back("agg");
    489                min.push_back(true);
    490                nRows.push_back(489);
    491                nCols.push_back(163);
    492                objValueTol.push_back(1.e-8);
    493                objValue.push_back(-3.5991767287e+07);
    494                bestStrategy.push_back(3);
    495                mpsName.push_back("agg2");
    496                min.push_back(true);
    497                nRows.push_back(517);
    498                nCols.push_back(302);
    499                objValueTol.push_back(1.e-8);
    500                objValue.push_back(-2.0239252356e+07);
    501                bestStrategy.push_back(3);
    502                mpsName.push_back("agg3");
    503                min.push_back(true);
    504                nRows.push_back(517);
    505                nCols.push_back(302);
    506                objValueTol.push_back(1.e-8);
    507                objValue.push_back(1.0312115935e+07);
    508                bestStrategy.push_back(4);
    509                mpsName.push_back("bandm");
    510                min.push_back(true);
    511                nRows.push_back(306);
    512                nCols.push_back(472);
    513                objValueTol.push_back(1.e-8);
    514                objValue.push_back(-1.5862801845e+02);
    515                bestStrategy.push_back(2);
    516                mpsName.push_back("beaconfd");
    517                min.push_back(true);
    518                nRows.push_back(174);
    519                nCols.push_back(262);
    520                objValueTol.push_back(1.e-8);
    521                objValue.push_back(3.3592485807e+04);
    522                bestStrategy.push_back(0);
    523                mpsName.push_back("bnl1");
    524                min.push_back(true);
    525                nRows.push_back(644);
    526                nCols.push_back(1175);
    527                objValueTol.push_back(1.e-8);
    528                objValue.push_back(1.9776295615E+03);
    529                bestStrategy.push_back(3);
    530                mpsName.push_back("bnl2");
    531                min.push_back(true);
    532                nRows.push_back(2325);
    533                nCols.push_back(3489);
    534                objValueTol.push_back(1.e-8);
    535                objValue.push_back(1.8112365404e+03);
    536                bestStrategy.push_back(3);
    537                mpsName.push_back("boeing1");
    538                min.push_back(true);
    539                nRows.push_back(/*351*/352);
    540                nCols.push_back(384);
    541                objValueTol.push_back(1.e-8);
    542                objValue.push_back(-3.3521356751e+02);
    543                bestStrategy.push_back(3);
    544                mpsName.push_back("boeing2");
    545                min.push_back(true);
    546                nRows.push_back(167);
    547                nCols.push_back(143);
    548                objValueTol.push_back(1.e-8);
    549                objValue.push_back(-3.1501872802e+02);
    550                bestStrategy.push_back(3);
    551                mpsName.push_back("bore3d");
    552                min.push_back(true);
    553                nRows.push_back(234);
    554                nCols.push_back(315);
    555                objValueTol.push_back(1.e-8);
    556                objValue.push_back(1.3730803942e+03);
    557                bestStrategy.push_back(3);
    558                mpsName.push_back("brandy");
    559                min.push_back(true);
    560                nRows.push_back(221);
    561                nCols.push_back(249);
    562                objValueTol.push_back(1.e-8);
    563                objValue.push_back(1.5185098965e+03);
    564                bestStrategy.push_back(3);
    565                mpsName.push_back("capri");
    566                min.push_back(true);
    567                nRows.push_back(272);
    568                nCols.push_back(353);
    569                objValueTol.push_back(1.e-8);
    570                objValue.push_back(2.6900129138e+03);
    571                bestStrategy.push_back(3);
    572                mpsName.push_back("cycle");
    573                min.push_back(true);
    574                nRows.push_back(1904);
    575                nCols.push_back(2857);
    576                objValueTol.push_back(1.e-9);
    577                objValue.push_back(-5.2263930249e+00);
    578                bestStrategy.push_back(3);
    579                mpsName.push_back("czprob");
    580                min.push_back(true);
    581                nRows.push_back(930);
    582                nCols.push_back(3523);
    583                objValueTol.push_back(1.e-8);
    584                objValue.push_back(2.1851966989e+06);
    585                bestStrategy.push_back(3);
    586                mpsName.push_back("d2q06c");
    587                min.push_back(true);
    588                nRows.push_back(2172);
    589                nCols.push_back(5167);
    590                objValueTol.push_back(1.e-7);
    591                objValue.push_back(122784.21557456);
    592                bestStrategy.push_back(0);
    593                mpsName.push_back("d6cube");
    594                min.push_back(true);
    595                nRows.push_back(416);
    596                nCols.push_back(6184);
    597                objValueTol.push_back(1.e-7);
    598                objValue.push_back(3.1549166667e+02);
    599                bestStrategy.push_back(3);
    600                mpsName.push_back("degen2");
    601                min.push_back(true);
    602                nRows.push_back(445);
    603                nCols.push_back(534);
    604                objValueTol.push_back(1.e-8);
    605                objValue.push_back(-1.4351780000e+03);
    606                bestStrategy.push_back(3);
    607                mpsName.push_back("degen3");
    608                min.push_back(true);
    609                nRows.push_back(1504);
    610                nCols.push_back(1818);
    611                objValueTol.push_back(1.e-8);
    612                objValue.push_back(-9.8729400000e+02);
    613                bestStrategy.push_back(2);
    614                mpsName.push_back("dfl001");
    615                min.push_back(true);
    616                nRows.push_back(6072);
    617                nCols.push_back(12230);
    618                objValueTol.push_back(1.e-5);
    619                objValue.push_back(1.1266396047E+07);
    620                bestStrategy.push_back(5);
    621                mpsName.push_back("e226");
    622                min.push_back(true);
    623                nRows.push_back(224);
    624                nCols.push_back(282);
    625                objValueTol.push_back(1.e-8);
    626                objValue.push_back(-1.8751929066e+01 + 7.113);
    627                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.
    628                mpsName.push_back("etamacro");
    629                min.push_back(true);
    630                nRows.push_back(401);
    631                nCols.push_back(688);
    632                objValueTol.push_back(1.e-6);
    633                objValue.push_back(-7.5571521774e+02 );
    634                bestStrategy.push_back(3);
    635                mpsName.push_back("fffff800");
    636                min.push_back(true);
    637                nRows.push_back(525);
    638                nCols.push_back(854);
    639                objValueTol.push_back(1.e-6);
    640                objValue.push_back(5.5567961165e+05);
    641                bestStrategy.push_back(3);
    642                mpsName.push_back("finnis");
    643                min.push_back(true);
    644                nRows.push_back(498);
    645                nCols.push_back(614);
    646                objValueTol.push_back(1.e-6);
    647                objValue.push_back(1.7279096547e+05);
    648                bestStrategy.push_back(3);
    649                mpsName.push_back("fit1d");
    650                min.push_back(true);
    651                nRows.push_back(25);
    652                nCols.push_back(1026);
    653                objValueTol.push_back(1.e-8);
    654                objValue.push_back(-9.1463780924e+03);
    655                bestStrategy.push_back(3 + 100);
    656                mpsName.push_back("fit1p");
    657                min.push_back(true);
    658                nRows.push_back(628);
    659                nCols.push_back(1677);
    660                objValueTol.push_back(1.e-8);
    661                objValue.push_back(9.1463780924e+03);
    662                bestStrategy.push_back(5 + 100);
    663                mpsName.push_back("fit2d");
    664                min.push_back(true);
    665                nRows.push_back(26);
    666                nCols.push_back(10500);
    667                objValueTol.push_back(1.e-8);
    668                objValue.push_back(-6.8464293294e+04);
    669                bestStrategy.push_back(3 + 100);
    670                mpsName.push_back("fit2p");
    671                min.push_back(true);
    672                nRows.push_back(3001);
    673                nCols.push_back(13525);
    674                objValueTol.push_back(1.e-9);
    675                objValue.push_back(6.8464293232e+04);
    676                bestStrategy.push_back(5 + 100);
    677                mpsName.push_back("forplan");
    678                min.push_back(true);
    679                nRows.push_back(162);
    680                nCols.push_back(421);
    681                objValueTol.push_back(1.e-6);
    682                objValue.push_back(-6.6421873953e+02);
    683                bestStrategy.push_back(3);
    684                mpsName.push_back("ganges");
    685                min.push_back(true);
    686                nRows.push_back(1310);
    687                nCols.push_back(1681);
    688                objValueTol.push_back(1.e-5);
    689                objValue.push_back(-1.0958636356e+05);
    690                bestStrategy.push_back(3);
    691                mpsName.push_back("gfrd-pnc");
    692                min.push_back(true);
    693                nRows.push_back(617);
    694                nCols.push_back(1092);
    695                objValueTol.push_back(1.e-8);
    696                objValue.push_back(6.9022359995e+06);
    697                bestStrategy.push_back(3);
    698                mpsName.push_back("greenbea");
    699                min.push_back(true);
    700                nRows.push_back(2393);
    701                nCols.push_back(5405);
    702                objValueTol.push_back(1.e-8);
    703                objValue.push_back(/*-7.2462405908e+07*/ -72555248.129846);
    704                bestStrategy.push_back(3);
    705                mpsName.push_back("greenbeb");
    706                min.push_back(true);
    707                nRows.push_back(2393);
    708                nCols.push_back(5405);
    709                objValueTol.push_back(1.e-8);
    710                objValue.push_back(/*-4.3021476065e+06*/ -4302260.2612066);
    711                bestStrategy.push_back(3);
    712                mpsName.push_back("grow15");
    713                min.push_back(true);
    714                nRows.push_back(301);
    715                nCols.push_back(645);
    716                objValueTol.push_back(1.e-8);
    717                objValue.push_back(-1.0687094129e+08);
    718                bestStrategy.push_back(4 + 100);
    719                mpsName.push_back("grow22");
    720                min.push_back(true);
    721                nRows.push_back(441);
    722                nCols.push_back(946);
    723                objValueTol.push_back(1.e-8);
    724                objValue.push_back(-1.6083433648e+08);
    725                bestStrategy.push_back(4 + 100);
    726                mpsName.push_back("grow7");
    727                min.push_back(true);
    728                nRows.push_back(141);
    729                nCols.push_back(301);
    730                objValueTol.push_back(1.e-8);
    731                objValue.push_back(-4.7787811815e+07);
    732                bestStrategy.push_back(4 + 100);
    733                mpsName.push_back("israel");
    734                min.push_back(true);
    735                nRows.push_back(175);
    736                nCols.push_back(142);
    737                objValueTol.push_back(1.e-8);
    738                objValue.push_back(-8.9664482186e+05);
    739                bestStrategy.push_back(2);
    740                mpsName.push_back("kb2");
    741                min.push_back(true);
    742                nRows.push_back(44);
    743                nCols.push_back(41);
    744                objValueTol.push_back(1.e-8);
    745                objValue.push_back(-1.7499001299e+03);
    746                bestStrategy.push_back(3);
    747                mpsName.push_back("lotfi");
    748                min.push_back(true);
    749                nRows.push_back(154);
    750                nCols.push_back(308);
    751                objValueTol.push_back(1.e-8);
    752                objValue.push_back(-2.5264706062e+01);
    753                bestStrategy.push_back(3);
    754                mpsName.push_back("maros");
    755                min.push_back(true);
    756                nRows.push_back(847);
    757                nCols.push_back(1443);
    758                objValueTol.push_back(1.e-8);
    759                objValue.push_back(-5.8063743701e+04);
    760                bestStrategy.push_back(3);
    761                mpsName.push_back("modszk1");
    762                min.push_back(true);
    763                nRows.push_back(688);
    764                nCols.push_back(1620);
    765                objValueTol.push_back(1.e-8);
    766                objValue.push_back(3.2061972906e+02);
    767                bestStrategy.push_back(3);
    768                mpsName.push_back("nesm");
    769                min.push_back(true);
    770                nRows.push_back(663);
    771                nCols.push_back(2923);
    772                objValueTol.push_back(1.e-5);
    773                objValue.push_back(1.4076073035e+07);
    774                bestStrategy.push_back(2);
    775                mpsName.push_back("perold");
    776                min.push_back(true);
    777                nRows.push_back(626);
    778                nCols.push_back(1376);
    779                objValueTol.push_back(1.e-6);
    780                objValue.push_back(-9.3807580773e+03);
    781                bestStrategy.push_back(3);
    782                //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);
    783                //mpsName.push_back("qap15");min.push_back(true);nRows.push_back(6331);nCols.push_back(22275);objValueTol.push_back(1.e-8);objValue.push_back(1.0409940410e+03);bestStrategy.push_back(3);
    784                mpsName.push_back("recipe");
    785                min.push_back(true);
    786                nRows.push_back(92);
    787                nCols.push_back(180);
    788                objValueTol.push_back(1.e-8);
    789                objValue.push_back(-2.6661600000e+02);
    790                bestStrategy.push_back(3);
    791                mpsName.push_back("sc105");
    792                min.push_back(true);
    793                nRows.push_back(106);
    794                nCols.push_back(103);
    795                objValueTol.push_back(1.e-8);
    796                objValue.push_back(-5.2202061212e+01);
    797                bestStrategy.push_back(3);
    798                mpsName.push_back("sc205");
    799                min.push_back(true);
    800                nRows.push_back(206);
    801                nCols.push_back(203);
    802                objValueTol.push_back(1.e-8);
    803                objValue.push_back(-5.2202061212e+01);
    804                bestStrategy.push_back(3);
    805                mpsName.push_back("sc50a");
    806                min.push_back(true);
    807                nRows.push_back(51);
    808                nCols.push_back(48);
    809                objValueTol.push_back(1.e-8);
    810                objValue.push_back(-6.4575077059e+01);
    811                bestStrategy.push_back(3);
    812                mpsName.push_back("sc50b");
    813                min.push_back(true);
    814                nRows.push_back(51);
    815                nCols.push_back(48);
    816                objValueTol.push_back(1.e-8);
    817                objValue.push_back(-7.0000000000e+01);
    818                bestStrategy.push_back(3);
    819                mpsName.push_back("scagr25");
    820                min.push_back(true);
    821                nRows.push_back(472);
    822                nCols.push_back(500);
    823                objValueTol.push_back(1.e-8);
    824                objValue.push_back(-1.4753433061e+07);
    825                bestStrategy.push_back(3);
    826                mpsName.push_back("scagr7");
    827                min.push_back(true);
    828                nRows.push_back(130);
    829                nCols.push_back(140);
    830                objValueTol.push_back(1.e-6);
    831                objValue.push_back(-2.3313892548e+06);
    832                bestStrategy.push_back(3);
    833                mpsName.push_back("scfxm1");
    834                min.push_back(true);
    835                nRows.push_back(331);
    836                nCols.push_back(457);
    837                objValueTol.push_back(1.e-8);
    838                objValue.push_back(1.8416759028e+04);
    839                bestStrategy.push_back(3);
    840                mpsName.push_back("scfxm2");
    841                min.push_back(true);
    842                nRows.push_back(661);
    843                nCols.push_back(914);
    844                objValueTol.push_back(1.e-8);
    845                objValue.push_back(3.6660261565e+04);
    846                bestStrategy.push_back(3);
    847                mpsName.push_back("scfxm3");
    848                min.push_back(true);
    849                nRows.push_back(991);
    850                nCols.push_back(1371);
    851                objValueTol.push_back(1.e-8);
    852                objValue.push_back(5.4901254550e+04);
    853                bestStrategy.push_back(3);
    854                mpsName.push_back("scorpion");
    855                min.push_back(true);
    856                nRows.push_back(389);
    857                nCols.push_back(358);
    858                objValueTol.push_back(1.e-8);
    859                objValue.push_back(1.8781248227e+03);
    860                bestStrategy.push_back(3);
    861                mpsName.push_back("scrs8");
    862                min.push_back(true);
    863                nRows.push_back(491);
    864                nCols.push_back(1169);
    865                objValueTol.push_back(1.e-5);
    866                objValue.push_back(9.0429998619e+02);
    867                bestStrategy.push_back(2);
    868                mpsName.push_back("scsd1");
    869                min.push_back(true);
    870                nRows.push_back(78);
    871                nCols.push_back(760);
    872                objValueTol.push_back(1.e-8);
    873                objValue.push_back(8.6666666743e+00);
    874                bestStrategy.push_back(3 + 100);
    875                mpsName.push_back("scsd6");
    876                min.push_back(true);
    877                nRows.push_back(148);
    878                nCols.push_back(1350);
    879                objValueTol.push_back(1.e-8);
    880                objValue.push_back(5.0500000078e+01);
    881                bestStrategy.push_back(3 + 100);
    882                mpsName.push_back("scsd8");
    883                min.push_back(true);
    884                nRows.push_back(398);
    885                nCols.push_back(2750);
    886                objValueTol.push_back(1.e-7);
    887                objValue.push_back(9.0499999993e+02);
    888                bestStrategy.push_back(1 + 100);
    889                mpsName.push_back("sctap1");
    890                min.push_back(true);
    891                nRows.push_back(301);
    892                nCols.push_back(480);
    893                objValueTol.push_back(1.e-8);
    894                objValue.push_back(1.4122500000e+03);
    895                bestStrategy.push_back(3);
    896                mpsName.push_back("sctap2");
    897                min.push_back(true);
    898                nRows.push_back(1091);
    899                nCols.push_back(1880);
    900                objValueTol.push_back(1.e-8);
    901                objValue.push_back(1.7248071429e+03);
    902                bestStrategy.push_back(3);
    903                mpsName.push_back("sctap3");
    904                min.push_back(true);
    905                nRows.push_back(1481);
    906                nCols.push_back(2480);
    907                objValueTol.push_back(1.e-8);
    908                objValue.push_back(1.4240000000e+03);
    909                bestStrategy.push_back(3);
    910                mpsName.push_back("seba");
    911                min.push_back(true);
    912                nRows.push_back(516);
    913                nCols.push_back(1028);
    914                objValueTol.push_back(1.e-8);
    915                objValue.push_back(1.5711600000e+04);
    916                bestStrategy.push_back(3);
    917                mpsName.push_back("share1b");
    918                min.push_back(true);
    919                nRows.push_back(118);
    920                nCols.push_back(225);
    921                objValueTol.push_back(1.e-8);
    922                objValue.push_back(-7.6589318579e+04);
    923                bestStrategy.push_back(3);
    924                mpsName.push_back("share2b");
    925                min.push_back(true);
    926                nRows.push_back(97);
    927                nCols.push_back(79);
    928                objValueTol.push_back(1.e-8);
    929                objValue.push_back(-4.1573224074e+02);
    930                bestStrategy.push_back(3);
    931                mpsName.push_back("shell");
    932                min.push_back(true);
    933                nRows.push_back(537);
    934                nCols.push_back(1775);
    935                objValueTol.push_back(1.e-8);
    936                objValue.push_back(1.2088253460e+09);
    937                bestStrategy.push_back(3);
    938                mpsName.push_back("ship04l");
    939                min.push_back(true);
    940                nRows.push_back(403);
    941                nCols.push_back(2118);
    942                objValueTol.push_back(1.e-8);
    943                objValue.push_back(1.7933245380e+06);
    944                bestStrategy.push_back(3);
    945                mpsName.push_back("ship04s");
    946                min.push_back(true);
    947                nRows.push_back(403);
    948                nCols.push_back(1458);
    949                objValueTol.push_back(1.e-8);
    950                objValue.push_back(1.7987147004e+06);
    951                bestStrategy.push_back(3);
    952                mpsName.push_back("ship08l");
    953                min.push_back(true);
    954                nRows.push_back(779);
    955                nCols.push_back(4283);
    956                objValueTol.push_back(1.e-8);
    957                objValue.push_back(1.9090552114e+06);
    958                bestStrategy.push_back(3);
    959                mpsName.push_back("ship08s");
    960                min.push_back(true);
    961                nRows.push_back(779);
    962                nCols.push_back(2387);
    963                objValueTol.push_back(1.e-8);
    964                objValue.push_back(1.9200982105e+06);
    965                bestStrategy.push_back(2);
    966                mpsName.push_back("ship12l");
    967                min.push_back(true);
    968                nRows.push_back(1152);
    969                nCols.push_back(5427);
    970                objValueTol.push_back(1.e-8);
    971                objValue.push_back(1.4701879193e+06);
    972                bestStrategy.push_back(3);
    973                mpsName.push_back("ship12s");
    974                min.push_back(true);
    975                nRows.push_back(1152);
    976                nCols.push_back(2763);
    977                objValueTol.push_back(1.e-8);
    978                objValue.push_back(1.4892361344e+06);
    979                bestStrategy.push_back(2);
    980                mpsName.push_back("sierra");
    981                min.push_back(true);
    982                nRows.push_back(1228);
    983                nCols.push_back(2036);
    984                objValueTol.push_back(1.e-8);
    985                objValue.push_back(1.5394362184e+07);
    986                bestStrategy.push_back(3);
    987                mpsName.push_back("stair");
    988                min.push_back(true);
    989                nRows.push_back(357);
    990                nCols.push_back(467);
    991                objValueTol.push_back(1.e-8);
    992                objValue.push_back(-2.5126695119e+02);
    993                bestStrategy.push_back(3);
    994                mpsName.push_back("standata");
    995                min.push_back(true);
    996                nRows.push_back(360);
    997                nCols.push_back(1075);
    998                objValueTol.push_back(1.e-8);
    999                objValue.push_back(1.2576995000e+03);
    1000                bestStrategy.push_back(3);
    1001                //mpsName.push_back("standgub");min.push_back(true);nRows.push_back(362);nCols.push_back(1184);objValueTol.push_back(1.e-8);objValue.push_back(1257.6995); bestStrategy.push_back(3);
    1002                mpsName.push_back("standmps");
    1003                min.push_back(true);
    1004                nRows.push_back(468);
    1005                nCols.push_back(1075);
    1006                objValueTol.push_back(1.e-8);
    1007                objValue.push_back(1.4060175000E+03);
    1008                bestStrategy.push_back(3);
    1009                mpsName.push_back("stocfor1");
    1010                min.push_back(true);
    1011                nRows.push_back(118);
    1012                nCols.push_back(111);
    1013                objValueTol.push_back(1.e-8);
    1014                objValue.push_back(-4.1131976219E+04);
    1015                bestStrategy.push_back(3);
    1016                mpsName.push_back("stocfor2");
    1017                min.push_back(true);
    1018                nRows.push_back(2158);
    1019                nCols.push_back(2031);
    1020                objValueTol.push_back(1.e-8);
    1021                objValue.push_back(-3.9024408538e+04);
    1022                bestStrategy.push_back(3);
    1023                //mpsName.push_back("stocfor3");min.push_back(true);nRows.push_back(16676);nCols.push_back(15695);objValueTol.push_back(1.e-8);objValue.push_back(-3.9976661576e+04);bestStrategy.push_back(3);
    1024                //mpsName.push_back("truss");min.push_back(true);nRows.push_back(1001);nCols.push_back(8806);objValueTol.push_back(1.e-8);objValue.push_back(4.5881584719e+05);bestStrategy.push_back(3);
    1025                mpsName.push_back("tuff");
    1026                min.push_back(true);
    1027                nRows.push_back(334);
    1028                nCols.push_back(587);
    1029                objValueTol.push_back(1.e-8);
    1030                objValue.push_back(2.9214776509e-01);
    1031                bestStrategy.push_back(3);
    1032                mpsName.push_back("vtpbase");
    1033                min.push_back(true);
    1034                nRows.push_back(199);
    1035                nCols.push_back(203);
    1036                objValueTol.push_back(1.e-8);
    1037                objValue.push_back(1.2983146246e+05);
    1038                bestStrategy.push_back(3);
    1039                mpsName.push_back("wood1p");
    1040                min.push_back(true);
    1041                nRows.push_back(245);
    1042                nCols.push_back(2594);
    1043                objValueTol.push_back(5.e-5);
    1044                objValue.push_back(1.4429024116e+00);
    1045                bestStrategy.push_back(3);
    1046                mpsName.push_back("woodw");
    1047                min.push_back(true);
    1048                nRows.push_back(1099);
    1049                nCols.push_back(8405);
    1050                objValueTol.push_back(1.e-8);
    1051                objValue.push_back(1.3044763331E+00);
    1052                bestStrategy.push_back(3);
    1053           } else {
    1054                // Just testing one
    1055                mpsName.push_back(empty.problemName());
    1056                min.push_back(true);
    1057                nRows.push_back(-1);
    1058                nCols.push_back(-1);
    1059                objValueTol.push_back(1.e-8);
    1060                objValue.push_back(0.0);
    1061                bestStrategy.push_back(0);
    1062                int iTest;
    1063                std::string alg;
    1064                for (iTest = 0; iTest < NUMBER_ALGORITHMS; iTest++) {
    1065                     ClpSolve solveOptions = setupForSolve(iTest, alg, 0);
    1066                     printf("%d %s ", iTest, alg.c_str());
    1067                     if (switchOff[iTest])
    1068                          printf("skipped by user\n");
    1069                     else if(solveOptions.getSolveType() == ClpSolve::notImplemented)
    1070                          printf("skipped as not available\n");
    1071                     else
    1072                          printf("will be tested\n");
    1073                }
    1074           }
    1075 
    1076           double timeTaken = 0.0;
    1077           if( !barrierAvailable)
    1078                switchOff[0] = 1;
    1079           // Loop once for each Mps File
    1080           for (m = 0; m < mpsName.size(); m++ ) {
    1081 #if FACTORIZATION_STATISTICS
    1082                if (nRows[m]<loSizeX||nRows[m]>=hiSizeX) {
    1083                  std::cerr << "  skipping mps file: " << mpsName[m]
    1084                            <<" as "<<nRows[m]
    1085                          << " (" << m + 1 << " out of " << mpsName.size() << ")" << std::endl;
    1086                  continue;
    1087                }
    1088 #endif
    1089                std::cerr << "  processing mps file: " << mpsName[m]
    1090                          << " (" << m + 1 << " out of " << mpsName.size() << ")" << std::endl;
    1091                AnySimplex solutionBase = empty;
    1092                std::string fn = dirNetlib + mpsName[m];
    1093                if (!empty.numberRows() || algorithm < 6) {
    1094                     // Read data mps file,
    1095                     CoinMpsIO mps;
    1096                     int nerrors = mps.readMps(fn.c_str(), "mps");
    1097                     if (nerrors) {
    1098                          std::cerr << "Error " << nerrors << " when reading model from "
    1099                                    << fn.c_str() << "! Aborting tests.\n";
    1100                          return 1;
    1101                     }
    1102                     solutionBase.loadProblem(*mps.getMatrixByCol(), mps.getColLower(),
    1103                                              mps.getColUpper(),
    1104                                              mps.getObjCoefficients(),
    1105                                              mps.getRowLower(), mps.getRowUpper());
    1106 
    1107                     solutionBase.setDblParam(ClpObjOffset, mps.objectiveOffset());
    1108                }
    1109 
    1110                // Runs through strategies
    1111                if (algorithm == 6 || algorithm == 7) {
    1112                     // algorithms tested are at top of file
    1113                     double testTime[NUMBER_ALGORITHMS];
    1114                     std::string alg[NUMBER_ALGORITHMS];
    1115                     int iTest;
    1116                     for (iTest = 0; iTest < NUMBER_ALGORITHMS; iTest++) {
    1117                          ClpSolve solveOptions = setupForSolve(iTest, alg[iTest], 1);
    1118                          if (solveOptions.getSolveType() != ClpSolve::notImplemented) {
    1119                               double time1 = CoinCpuTime();
    1120                               ClpSimplex solution = solutionBase;
    1121                               if (solution.maximumSeconds() < 0.0)
    1122                                    solution.setMaximumSeconds(120.0);
    1123                               if (doVector) {
    1124                                    ClpMatrixBase * matrix = solution.clpMatrix();
    1125                                    if (dynamic_cast< ClpPackedMatrix*>(matrix)) {
    1126                                         ClpPackedMatrix * clpMatrix = dynamic_cast< ClpPackedMatrix*>(matrix);
    1127                                         clpMatrix->makeSpecialColumnCopy();
    1128                                    }
    1129                               }
    1130                               solution.initialSolve(solveOptions);
    1131                               double time2 = CoinCpuTime() - time1;
    1132                               testTime[iTest] = time2;
    1133                               printf("Finished %s Took %g seconds (%d iterations) - status %d\n",
    1134                                      mpsName[m].c_str(),time2, solution.problemStatus(),solution.numberIterations());
    1135                               if (solution.problemStatus())
    1136                                    testTime[iTest] = 1.0e20;
    1137                          } else {
    1138                               testTime[iTest] = 1.0e30;
    1139                          }
    1140                     }
    1141                     int iBest = -1;
    1142                     double dBest = 1.0e10;
    1143                     printf("%s", fn.c_str());
    1144                     for (iTest = 0; iTest < NUMBER_ALGORITHMS; iTest++) {
    1145                          if (testTime[iTest] < 1.0e30) {
    1146                               printf(" %s %g",
    1147                                      alg[iTest].c_str(), testTime[iTest]);
    1148                               if (testTime[iTest] < dBest) {
    1149                                    dBest = testTime[iTest];
    1150                                    iBest = iTest;
    1151                               }
    1152                          }
    1153                     }
    1154                     printf("\n");
    1155                     if (iBest >= 0)
    1156                          printf("Best strategy for %s is %s (%d) which takes %g seconds\n",
    1157                                 fn.c_str(), alg[iBest].c_str(), iBest, testTime[iBest]);
    1158                     else
    1159                          printf("No strategy finished in time limit\n");
    1160                     continue;
    1161                }
    1162                double time1 = CoinCpuTime();
    1163                AnySimplex solution = solutionBase;
     1153        }
     1154        printf("\n");
     1155        if (iBest >= 0)
     1156          printf("Best strategy for %s is %s (%d) which takes %g seconds\n",
     1157            fn.c_str(), alg[iBest].c_str(), iBest, testTime[iBest]);
     1158        else
     1159          printf("No strategy finished in time limit\n");
     1160        continue;
     1161      }
     1162      double time1 = CoinCpuTime();
     1163      AnySimplex solution = solutionBase;
    11641164#if 0
    11651165               solution.setOptimizationDirection(-1);
     
    11721172               }
    11731173#endif
    1174                ClpSolve::SolveType method;
    1175                ClpSolve solveOptions = solveOptionsIn;
    1176                std::string nameAlgorithm;
    1177                if (algorithm != 5) {
    1178                     if (algorithm == 0) {
    1179                          method = ClpSolve::useDual;
    1180                          nameAlgorithm = "dual";
    1181                     } else if (algorithm == 1) {
    1182                          method = ClpSolve::usePrimalorSprint;
    1183                          nameAlgorithm = "primal";
    1184                     } else if (algorithm == 3) {
    1185                          method = ClpSolve::automatic;
    1186                          nameAlgorithm = "either";
    1187                     } else {
    1188                          nameAlgorithm = "barrier-slow";
     1174      ClpSolve::SolveType method;
     1175      ClpSolve solveOptions = solveOptionsIn;
     1176      std::string nameAlgorithm;
     1177      if (algorithm != 5) {
     1178        if (algorithm == 0) {
     1179          method = ClpSolve::useDual;
     1180          nameAlgorithm = "dual";
     1181        } else if (algorithm == 1) {
     1182          method = ClpSolve::usePrimalorSprint;
     1183          nameAlgorithm = "primal";
     1184        } else if (algorithm == 3) {
     1185          method = ClpSolve::automatic;
     1186          nameAlgorithm = "either";
     1187        } else {
     1188          nameAlgorithm = "barrier-slow";
    11891189#if defined(COIN_HAS_AMD) || defined(COIN_HAS_CHOLMOD) || defined(COIN_HAS_GLPK)
    1190                          solveOptions.setSpecialOption(4, 4);
    1191                          nameAlgorithm = "barrier-UFL";
     1190          solveOptions.setSpecialOption(4, 4);
     1191          nameAlgorithm = "barrier-UFL";
    11921192#endif
    11931193#ifdef COIN_HAS_WSMP
    1194                          solveOptions.setSpecialOption(4, 2);
    1195                          nameAlgorithm = "barrier-WSSMP";
     1194          solveOptions.setSpecialOption(4, 2);
     1195          nameAlgorithm = "barrier-WSSMP";
    11961196#endif
    11971197#ifdef COIN_HAS_MUMPS
    1198                          solveOptions.setSpecialOption(4, 6);
    1199                          nameAlgorithm = "barrier-MUMPS";
    1200 #endif
    1201                          method = ClpSolve::useBarrier;
    1202                     }
    1203                     solveOptions.setSolveType(method);
    1204                } else {
    1205                     int iAlg = bestStrategy[m];
    1206                     int presolveOff = iAlg / 100;
    1207                     iAlg = iAlg % 100;
    1208                     if( !barrierAvailable && iAlg == 0) {
    1209                          if (nRows[m] != 2172)
    1210                               iAlg = 5; // try primal
    1211                          else
    1212                               iAlg = 3; // d2q06c
    1213                     }
    1214                     solveOptions = setupForSolve(iAlg, nameAlgorithm, 0);
    1215                     if (presolveOff)
    1216                          solveOptions.setPresolveType(ClpSolve::presolveOff);
    1217                }
    1218                if (doVector) {
    1219                     ClpMatrixBase * matrix = solution.clpMatrix();
    1220                     if (dynamic_cast< ClpPackedMatrix*>(matrix)) {
    1221                          ClpPackedMatrix * clpMatrix = dynamic_cast< ClpPackedMatrix*>(matrix);
    1222                          clpMatrix->makeSpecialColumnCopy();
    1223                     }
    1224                }
     1198          solveOptions.setSpecialOption(4, 6);
     1199          nameAlgorithm = "barrier-MUMPS";
     1200#endif
     1201          method = ClpSolve::useBarrier;
     1202        }
     1203        solveOptions.setSolveType(method);
     1204      } else {
     1205        int iAlg = bestStrategy[m];
     1206        int presolveOff = iAlg / 100;
     1207        iAlg = iAlg % 100;
     1208        if (!barrierAvailable && iAlg == 0) {
     1209          if (nRows[m] != 2172)
     1210            iAlg = 5; // try primal
     1211          else
     1212            iAlg = 3; // d2q06c
     1213        }
     1214        solveOptions = setupForSolve(iAlg, nameAlgorithm, 0);
     1215        if (presolveOff)
     1216          solveOptions.setPresolveType(ClpSolve::presolveOff);
     1217      }
     1218      if (doVector) {
     1219        ClpMatrixBase *matrix = solution.clpMatrix();
     1220        if (dynamic_cast< ClpPackedMatrix * >(matrix)) {
     1221          ClpPackedMatrix *clpMatrix = dynamic_cast< ClpPackedMatrix * >(matrix);
     1222          clpMatrix->makeSpecialColumnCopy();
     1223        }
     1224      }
    12251225#if FACTORIZATION_STATISTICS
    1226                double timesOne[3*3*3*3*3];
    1227                memset(timesOne,0,sizeof(timesOne));
    1228                int iterationsOne[3*3*3*3*3];
    1229                memset(iterationsOne,0,sizeof(iterationsOne));
    1230                AnySimplex saveModel(solution);
    1231                double time2;
     1226      double timesOne[3 * 3 * 3 * 3 * 3];
     1227      memset(timesOne, 0, sizeof(timesOne));
     1228      int iterationsOne[3 * 3 * 3 * 3 * 3];
     1229      memset(iterationsOne, 0, sizeof(iterationsOne));
     1230      AnySimplex saveModel(solution);
     1231      double time2;
    12321232#if 0
    12331233               solution.initialSolve(solveOptions);
     
    12461246#define hiD 2
    12471247#define hiE 2
    1248                time2 = CoinCpuTime();
    1249                for (int iA=loA;iA<hiA;iA++) {
    1250                  *externalA=useA[iA];
    1251                  for (int iB=loB;iB<hiB;iB++) {
    1252                    *externalB=useB[iB];
    1253                    for (int iC=loC;iC<hiC;iC++) {
    1254                      *externalC=useC[iC];
    1255                      for (int iD=loD;iD<hiD;iD++) {
    1256                        *externalD=useD[iD];
    1257                        for (int iE=loE;iE<hiE;iE++) {
    1258                          *externalE=useE[iE];
    1259                          solution=saveModel;
    1260                         solution.initialSolve(solveOptions);
    1261                          double time3=CoinCpuTime();
    1262                          printf("Finished %s Took %g seconds (%d iterations) - status %d\n",
    1263                                 mpsName[m].c_str(),time3-time2,solution.numberIterations(), solution.problemStatus());
    1264                          iterationsOne[iA+3*iB+9*iC+27*iD+81*iE]=solution.numberIterations();
    1265                          timesOne[iA+3*iB+9*iC+27*iD+81*iE]=time3-time2;
    1266                          times[iA+3*iB+9*iC+27*iD+81*iE]+=time3-time2;
    1267                          time2=time3;
    1268                        }
    1269                      }
    1270                    }
    1271                 }
    1272                }
    1273                double bestTime=1.0e100;
    1274                int iBestTime=-1;
    1275                double bestTimePer=1.0e100;
    1276                int iBestTimePer=-1;
    1277                int bestIterations=1000000;
    1278                int iBestIterations=-1;
    1279                printf("TTimes ");
    1280                for (int i=0;i<3*3*3*3*3;i++) {
    1281                 // skip if not done
    1282                 if (!iterationsOne[i])
    1283                    continue;
    1284                  double average=timesOne[i]/static_cast<double>(iterationsOne[i]);
    1285                  if (timesOne[i]<bestTime) {
    1286                    bestTime=timesOne[i];
    1287                    iBestTime=i;
    1288                 }
    1289                  if (average<bestTimePer) {
    1290                    bestTimePer=average;
    1291                    iBestTimePer=i;
    1292                 }
    1293                  if (iterationsOne[i]<bestIterations) {
    1294                    bestIterations=iterationsOne[i];
    1295                    iBestIterations=i;
    1296                 }
    1297                  printf("%.2f ",timesOne[i]);
    1298                }
    1299                printf("\n");
    1300                int iA,iB,iC,iD,iE;
    1301                iA=iBestTime;
    1302                iE=iA/81;
    1303                iA-=81*iE;
    1304                iD=iA/27;
    1305                iA-=27*iD;
    1306                iC=iA/9;
    1307                iA-=9*iC;
    1308                iB=iA/3;
    1309                iA-=3*iB;
    1310                printf("Best time %.2f %s=%g %s=%g %s=%g %s=%g %s=%g\n",bestTime,choiceA,useA[iA],
    1311                       choiceB,useB[iB],choiceC,useC[iC],choiceD,useD[iD],choiceE,useE[iE]);
    1312                iA=iBestTimePer;
    1313                iE=iA/81;
    1314                iA-=81*iE;
    1315                iD=iA/27;
    1316                iA-=27*iD;
    1317                iC=iA/9;
    1318                iA-=9*iC;
    1319                iB=iA/3;
    1320                iA-=3*iB;
    1321                printf("Best time per iteration %g %s=%g %s=%g %s=%g %s=%g %s=%g\n",bestTimePer,choiceA,useA[iA],
    1322                       choiceB,useB[iB],choiceC,useC[iC],choiceD,useD[iD],choiceE,useE[iE]);
    1323                iA=iBestIterations;
    1324                iE=iA/81;
    1325                iA-=81*iE;
    1326                iD=iA/27;
    1327                iA-=27*iD;
    1328                iC=iA/9;
    1329                iA-=9*iC;
    1330                iB=iA/3;
    1331                iA-=3*iB;
    1332                printf("Best iterations %d %s=%g %s=%g %s=%g %s=%g %s=%g\n",bestIterations,choiceA,useA[iA],
    1333                       choiceB,useB[iB],choiceC,useC[iC],choiceD,useD[iD],choiceE,useE[iE]);
     1248      time2 = CoinCpuTime();
     1249      for (int iA = loA; iA < hiA; iA++) {
     1250        *externalA = useA[iA];
     1251        for (int iB = loB; iB < hiB; iB++) {
     1252          *externalB = useB[iB];
     1253          for (int iC = loC; iC < hiC; iC++) {
     1254            *externalC = useC[iC];
     1255            for (int iD = loD; iD < hiD; iD++) {
     1256              *externalD = useD[iD];
     1257              for (int iE = loE; iE < hiE; iE++) {
     1258                *externalE = useE[iE];
     1259                solution = saveModel;
     1260                solution.initialSolve(solveOptions);
     1261                double time3 = CoinCpuTime();
     1262                printf("Finished %s Took %g seconds (%d iterations) - status %d\n",
     1263                  mpsName[m].c_str(), time3 - time2, solution.numberIterations(), solution.problemStatus());
     1264                iterationsOne[iA + 3 * iB + 9 * iC + 27 * iD + 81 * iE] = solution.numberIterations();
     1265                timesOne[iA + 3 * iB + 9 * iC + 27 * iD + 81 * iE] = time3 - time2;
     1266                times[iA + 3 * iB + 9 * iC + 27 * iD + 81 * iE] += time3 - time2;
     1267                time2 = time3;
     1268              }
     1269            }
     1270          }
     1271        }
     1272      }
     1273      double bestTime = 1.0e100;
     1274      int iBestTime = -1;
     1275      double bestTimePer = 1.0e100;
     1276      int iBestTimePer = -1;
     1277      int bestIterations = 1000000;
     1278      int iBestIterations = -1;
     1279      printf("TTimes ");
     1280      for (int i = 0; i < 3 * 3 * 3 * 3 * 3; i++) {
     1281        // skip if not done
     1282        if (!iterationsOne[i])
     1283          continue;
     1284        double average = timesOne[i] / static_cast< double >(iterationsOne[i]);
     1285        if (timesOne[i] < bestTime) {
     1286          bestTime = timesOne[i];
     1287          iBestTime = i;
     1288        }
     1289        if (average < bestTimePer) {
     1290          bestTimePer = average;
     1291          iBestTimePer = i;
     1292        }
     1293        if (iterationsOne[i] < bestIterations) {
     1294          bestIterations = iterationsOne[i];
     1295          iBestIterations = i;
     1296        }
     1297        printf("%.2f ", timesOne[i]);
     1298      }
     1299      printf("\n");
     1300      int iA, iB, iC, iD, iE;
     1301      iA = iBestTime;
     1302      iE = iA / 81;
     1303      iA -= 81 * iE;
     1304      iD = iA / 27;
     1305      iA -= 27 * iD;
     1306      iC = iA / 9;
     1307      iA -= 9 * iC;
     1308      iB = iA / 3;
     1309      iA -= 3 * iB;
     1310      printf("Best time %.2f %s=%g %s=%g %s=%g %s=%g %s=%g\n", bestTime, choiceA, useA[iA],
     1311        choiceB, useB[iB], choiceC, useC[iC], choiceD, useD[iD], choiceE, useE[iE]);
     1312      iA = iBestTimePer;
     1313      iE = iA / 81;
     1314      iA -= 81 * iE;
     1315      iD = iA / 27;
     1316      iA -= 27 * iD;
     1317      iC = iA / 9;
     1318      iA -= 9 * iC;
     1319      iB = iA / 3;
     1320      iA -= 3 * iB;
     1321      printf("Best time per iteration %g %s=%g %s=%g %s=%g %s=%g %s=%g\n", bestTimePer, choiceA, useA[iA],
     1322        choiceB, useB[iB], choiceC, useC[iC], choiceD, useD[iD], choiceE, useE[iE]);
     1323      iA = iBestIterations;
     1324      iE = iA / 81;
     1325      iA -= 81 * iE;
     1326      iD = iA / 27;
     1327      iA -= 27 * iD;
     1328      iC = iA / 9;
     1329      iA -= 9 * iC;
     1330      iB = iA / 3;
     1331      iA -= 3 * iB;
     1332      printf("Best iterations %d %s=%g %s=%g %s=%g %s=%g %s=%g\n", bestIterations, choiceA, useA[iA],
     1333        choiceB, useB[iB], choiceC, useC[iC], choiceD, useD[iD], choiceE, useE[iE]);
    13341334#else
    1335                solution.initialSolve(solveOptions);
    1336                double time2 = CoinCpuTime() - time1;
    1337                timeTaken += time2;
    1338                printf("%s took %g seconds using algorithm %s\n", fn.c_str(), time2, nameAlgorithm.c_str());
    1339 #endif
    1340                // Test objective solution value
    1341                {
    1342                     if (!solution.isProvenOptimal()) {
    1343                       std::cerr <<"** NOT OPTIMAL ";
    1344                       numberFailures++;
    1345                     }
    1346                     double soln = solution.objectiveValue();
    1347                     CoinRelFltEq eq(objValueTol[m]);
    1348                     std::cerr << soln << ",  " << objValue[m] << " diff " <<
    1349                               soln - objValue[m] << std::endl;
    1350                     if(!eq(soln, objValue[m])) {
    1351                          printf("** difference fails\n");
    1352                          numberFailures++;
    1353                     }
    1354                }
     1335      solution.initialSolve(solveOptions);
     1336      double time2 = CoinCpuTime() - time1;
     1337      timeTaken += time2;
     1338      printf("%s took %g seconds using algorithm %s\n", fn.c_str(), time2, nameAlgorithm.c_str());
     1339#endif
     1340      // Test objective solution value
     1341      {
     1342        if (!solution.isProvenOptimal()) {
     1343          std::cerr << "** NOT OPTIMAL ";
     1344          numberFailures++;
     1345        }
     1346        double soln = solution.objectiveValue();
     1347        CoinRelFltEq eq(objValueTol[m]);
     1348        std::cerr << soln << ",  " << objValue[m] << " diff " << soln - objValue[m] << std::endl;
     1349        if (!eq(soln, objValue[m])) {
     1350          printf("** difference fails\n");
     1351          numberFailures++;
     1352        }
     1353      }
     1354    }
     1355    printf("Total time %g seconds\n", timeTaken);
     1356#if FACTORIZATION_STATISTICS
     1357    double bestTime = 1.0e100;
     1358    int iBestTime = -1;
     1359    for (int i = 0; i < 3 * 3 * 3 * 3 * 3; i++) {
     1360      if (times[i] && times[i] < bestTime) {
     1361        bestTime = times[i];
     1362        iBestTime = i;
     1363      }
     1364    }
     1365    for (int iA = 0; iA < hiA; iA++) {
     1366      for (int iB = 0; iB < hiB; iB++) {
     1367        for (int iC = 0; iC < hiC; iC++) {
     1368          for (int iD = loD; iD < hiD; iD++) {
     1369            for (int iE = loE; iE < hiE; iE++) {
     1370              int k = iA + 3 * iB + 9 * iC + 27 * iD + 81 * iE;
     1371              double thisTime = times[k];
     1372              if (thisTime) {
     1373                if (k == iBestTime)
     1374                  printf("** ");
     1375                printf("%s=%g %s=%g %s=%g %s=%g %s=%g - time %.2f\n",
     1376                  choiceA, useA[iA],
     1377                  choiceB, useB[iB], choiceC, useC[iC],
     1378                  choiceD, useD[iD], choiceE, useE[iE], thisTime);
     1379              }
     1380            }
    13551381          }
    1356           printf("Total time %g seconds\n", timeTaken);
    1357 #if FACTORIZATION_STATISTICS
    1358           double bestTime=1.0e100;
    1359           int iBestTime=-1;
    1360           for (int i=0;i<3*3*3*3*3;i++) {
    1361             if (times[i]&&times[i]<bestTime) {
    1362               bestTime=times[i];
    1363               iBestTime=i;
    1364             }
    1365           }
    1366           for (int iA=0;iA<hiA;iA++) {
    1367             for (int iB=0;iB<hiB;iB++) {
    1368               for (int iC=0;iC<hiC;iC++) {
    1369                 for (int iD=loD;iD<hiD;iD++) {
    1370                   for (int iE=loE;iE<hiE;iE++) {
    1371                     int k=iA+3*iB+9*iC+27*iD+81*iE;
    1372                     double thisTime=times[k];
    1373                     if (thisTime) {
    1374                       if (k==iBestTime)
    1375                         printf("** ");
    1376                       printf("%s=%g %s=%g %s=%g %s=%g %s=%g - time %.2f\n",
    1377                              choiceA,useA[iA],
    1378                              choiceB,useB[iB],choiceC,useC[iC],
    1379                              choiceD,useD[iD],choiceE,useE[iE],thisTime);
    1380                     }
    1381                   }
    1382                 }
    1383               }
    1384             }
    1385           }
    1386           //exit(0);
    1387 #endif
    1388      } else {
    1389           testingMessage( "***Skipped Testing on netlib    ***\n" );
    1390           testingMessage( "***use -netlib to test class***\n" );
    1391      }
    1392      if (!numberFailures)
    1393        testingMessage( "All tests completed successfully\n" );
    1394      else
    1395        testingMessage( "Some tests failed\n" );
    1396      return 0;
     1382        }
     1383      }
     1384    }
     1385    //exit(0);
     1386#endif
     1387  } else {
     1388    testingMessage("***Skipped Testing on netlib    ***\n");
     1389    testingMessage("***use -netlib to test class***\n");
     1390  }
     1391  if (!numberFailures)
     1392    testingMessage("All tests completed successfully\n");
     1393  else
     1394    testingMessage("Some tests failed\n");
     1395  return 0;
    13971396}
    13981397
    1399 
    14001398// Display message on stdout and stderr
    1401 void testingMessage( const char * const msg )
     1399void testingMessage(const char *const msg)
    14021400{
    1403      std::cerr << msg;
    1404      //cout <<endl <<"*****************************************"
    1405      //     <<endl <<msg <<endl;
     1401  std::cerr << msg;
     1402  //cout <<endl <<"*****************************************"
     1403  //     <<endl <<msg <<endl;
    14061404}
    14071405
    14081406//--------------------------------------------------------------------------
    14091407// test factorization methods and simplex method and simple barrier
    1410 void
    1411 ClpSimplexUnitTest(const std::string & dirSample)
     1408void ClpSimplexUnitTest(const std::string &dirSample)
    14121409{
    14131410
    1414      CoinRelFltEq eq(0.000001);
    1415 
    1416      {
    1417           ClpSimplex solution;
    1418 
    1419           // matrix data
    1420           //deliberate hiccup of 2 between 0 and 1
    1421           CoinBigIndex start[5] = {0, 4, 7, 8, 9};
    1422           int length[5] = {2, 3, 1, 1, 1};
    1423           int rows[11] = {0, 2, -1, -1, 0, 1, 2, 0, 1, 2};
    1424           double elements[11] = {7.0, 2.0, 1.0e10, 1.0e10, -2.0, 1.0, -2.0, 1, 1, 1};
    1425           CoinPackedMatrix matrix(true, 3, 5, 8, elements, rows, start, length);
    1426 
    1427           // rim data
    1428           double objective[7] = { -4.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0};
    1429           double rowLower[5] = {14.0, 3.0, 3.0, 1.0e10, 1.0e10};
    1430           double rowUpper[5] = {14.0, 3.0, 3.0, -1.0e10, -1.0e10};
    1431           double colLower[7] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
    1432           double colUpper[7] = {100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0};
    1433 
    1434           // basis 1
    1435           int rowBasis1[3] = { -1, -1, -1};
    1436           int colBasis1[5] = {1, 1, -1, -1, 1};
    1437           solution.loadProblem(matrix, colLower, colUpper, objective,
    1438                                rowLower, rowUpper);
    1439           int i;
    1440           solution.createStatus();
    1441           for (i = 0; i < 3; i++) {
    1442                if (rowBasis1[i] < 0) {
    1443                     solution.setRowStatus(i, ClpSimplex::atLowerBound);
    1444                } else {
    1445                     solution.setRowStatus(i, ClpSimplex::basic);
    1446                }
    1447           }
    1448           for (i = 0; i < 5; i++) {
    1449                if (colBasis1[i] < 0) {
    1450                     solution.setColumnStatus(i, ClpSimplex::atLowerBound);
    1451                } else {
    1452                     solution.setColumnStatus(i, ClpSimplex::basic);
    1453                }
    1454           }
    1455           solution.setLogLevel(3 + 4 + 8 + 16 + 32);
    1456           solution.primal();
    1457           for (i = 0; i < 3; i++) {
    1458                if (rowBasis1[i] < 0) {
    1459                     solution.setRowStatus(i, ClpSimplex::atLowerBound);
    1460                } else {
    1461                     solution.setRowStatus(i, ClpSimplex::basic);
    1462                }
    1463           }
    1464           for (i = 0; i < 5; i++) {
    1465                if (colBasis1[i] < 0) {
    1466                     solution.setColumnStatus(i, ClpSimplex::atLowerBound);
    1467                } else {
    1468                     solution.setColumnStatus(i, ClpSimplex::basic);
    1469                }
    1470           }
    1471           // intricate stuff does not work with scaling
    1472           solution.scaling(0);
     1411  CoinRelFltEq eq(0.000001);
     1412
     1413  {
     1414    ClpSimplex solution;
     1415
     1416    // matrix data
     1417    //deliberate hiccup of 2 between 0 and 1
     1418    CoinBigIndex start[5] = { 0, 4, 7, 8, 9 };
     1419    int length[5] = { 2, 3, 1, 1, 1 };
     1420    int rows[11] = { 0, 2, -1, -1, 0, 1, 2, 0, 1, 2 };
     1421    double elements[11] = { 7.0, 2.0, 1.0e10, 1.0e10, -2.0, 1.0, -2.0, 1, 1, 1 };
     1422    CoinPackedMatrix matrix(true, 3, 5, 8, elements, rows, start, length);
     1423
     1424    // rim data
     1425    double objective[7] = { -4.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0 };
     1426    double rowLower[5] = { 14.0, 3.0, 3.0, 1.0e10, 1.0e10 };
     1427    double rowUpper[5] = { 14.0, 3.0, 3.0, -1.0e10, -1.0e10 };
     1428    double colLower[7] = { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 };
     1429    double colUpper[7] = { 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0 };
     1430
     1431    // basis 1
     1432    int rowBasis1[3] = { -1, -1, -1 };
     1433    int colBasis1[5] = { 1, 1, -1, -1, 1 };
     1434    solution.loadProblem(matrix, colLower, colUpper, objective,
     1435      rowLower, rowUpper);
     1436    int i;
     1437    solution.createStatus();
     1438    for (i = 0; i < 3; i++) {
     1439      if (rowBasis1[i] < 0) {
     1440        solution.setRowStatus(i, ClpSimplex::atLowerBound);
     1441      } else {
     1442        solution.setRowStatus(i, ClpSimplex::basic);
     1443      }
     1444    }
     1445    for (i = 0; i < 5; i++) {
     1446      if (colBasis1[i] < 0) {
     1447        solution.setColumnStatus(i, ClpSimplex::atLowerBound);
     1448      } else {
     1449        solution.setColumnStatus(i, ClpSimplex::basic);
     1450      }
     1451    }
     1452    solution.setLogLevel(3 + 4 + 8 + 16 + 32);
     1453    solution.primal();
     1454    for (i = 0; i < 3; i++) {
     1455      if (rowBasis1[i] < 0) {
     1456        solution.setRowStatus(i, ClpSimplex::atLowerBound);
     1457      } else {
     1458        solution.setRowStatus(i, ClpSimplex::basic);
     1459      }
     1460    }
     1461    for (i = 0; i < 5; i++) {
     1462      if (colBasis1[i] < 0) {
     1463        solution.setColumnStatus(i, ClpSimplex::atLowerBound);
     1464      } else {
     1465        solution.setColumnStatus(i, ClpSimplex::basic);
     1466      }
     1467    }
     1468    // intricate stuff does not work with scaling
     1469    solution.scaling(0);
    14731470#ifndef NDEBUG
    1474           int returnCode = solution.factorize ( );
    1475           assert(!returnCode);
     1471    int returnCode = solution.factorize();
     1472    assert(!returnCode);
    14761473#else
    1477           solution.factorize ( );
    1478 #endif
    1479           const double * colsol = solution.primalColumnSolution();
    1480           const double * rowsol = solution.primalRowSolution();
    1481           solution.getSolution(rowsol, colsol);
     1474    solution.factorize();
     1475#endif
     1476    const double *colsol = solution.primalColumnSolution();
     1477    const double *rowsol = solution.primalRowSolution();
     1478    solution.getSolution(rowsol, colsol);
    14821479#ifndef NDEBUG
    1483           double colsol1[5] = {20.0 / 7.0, 3.0, 0.0, 0.0, 23.0 / 7.0};
    1484           for (i = 0; i < 5; i++) {
    1485                assert(eq(colsol[i], colsol1[i]));
    1486           }
    1487           // now feed in again without actually doing factorization
    1488 #endif
    1489           ClpFactorization factorization2 = *solution.factorization();
    1490           ClpSimplex solution2 = solution;
    1491           solution2.setFactorization(factorization2);
    1492           solution2.createStatus();
    1493           for (i = 0; i < 3; i++) {
    1494                if (rowBasis1[i] < 0) {
    1495                     solution2.setRowStatus(i, ClpSimplex::atLowerBound);
    1496                } else {
    1497                     solution2.setRowStatus(i, ClpSimplex::basic);
    1498                }
    1499           }
    1500           for (i = 0; i < 5; i++) {
    1501                if (colBasis1[i] < 0) {
    1502                     solution2.setColumnStatus(i, ClpSimplex::atLowerBound);
    1503                } else {
    1504                     solution2.setColumnStatus(i, ClpSimplex::basic);
    1505                }
    1506           }
    1507           // intricate stuff does not work with scaling
    1508           solution2.scaling(0);
    1509           solution2.getSolution(rowsol, colsol);
    1510           colsol = solution2.primalColumnSolution();
    1511           rowsol = solution2.primalRowSolution();
    1512           for (i = 0; i < 5; i++) {
    1513                assert(eq(colsol[i], colsol1[i]));
    1514           }
    1515           solution2.setDualBound(0.1);
    1516           solution2.dual();
    1517           objective[2] = -1.0;
    1518           objective[3] = -0.5;
    1519           objective[4] = 10.0;
    1520           solution.dual();
    1521           for (i = 0; i < 3; i++) {
    1522                rowLower[i] = -1.0e20;
    1523                colUpper[i+2] = 0.0;
    1524           }
    1525           solution.setLogLevel(3);
    1526           solution.dual();
    1527           double rowObjective[] = {1.0, 0.5, -10.0};
    1528           solution.loadProblem(matrix, colLower, colUpper, objective,
    1529                                rowLower, rowUpper, rowObjective);
    1530           solution.dual();
    1531           solution.loadProblem(matrix, colLower, colUpper, objective,
    1532                                rowLower, rowUpper, rowObjective);
    1533           solution.primal();
    1534      }
     1480    double colsol1[5] = { 20.0 / 7.0, 3.0, 0.0, 0.0, 23.0 / 7.0 };
     1481    for (i = 0; i < 5; i++) {
     1482      assert(eq(colsol[i], colsol1[i]));
     1483    }
     1484    // now feed in again without actually doing factorization
     1485#endif
     1486    ClpFactorization factorization2 = *solution.factorization();
     1487    ClpSimplex solution2 = solution;
     1488    solution2.setFactorization(factorization2);
     1489    solution2.createStatus();
     1490    for (i = 0; i < 3; i++) {
     1491      if (rowBasis1[i] < 0) {
     1492        solution2.setRowStatus(i, ClpSimplex::atLowerBound);
     1493      } else {
     1494        solution2.setRowStatus(i, ClpSimplex::basic);
     1495      }
     1496    }
     1497    for (i = 0; i < 5; i++) {
     1498      if (colBasis1[i] < 0) {
     1499        solution2.setColumnStatus(i, ClpSimplex::atLowerBound);
     1500      } else {
     1501        solution2.setColumnStatus(i, ClpSimplex::basic);
     1502      }
     1503    }
     1504    // intricate stuff does not work with scaling
     1505    solution2.scaling(0);
     1506    solution2.getSolution(rowsol, colsol);
     1507    colsol = solution2.primalColumnSolution();
     1508    rowsol = solution2.primalRowSolution();
     1509    for (i = 0; i < 5; i++) {
     1510      assert(eq(colsol[i], colsol1[i]));
     1511    }
     1512    solution2.setDualBound(0.1);
     1513    solution2.dual();
     1514    objective[2] = -1.0;
     1515    objective[3] = -0.5;
     1516    objective[4] = 10.0;
     1517    solution.dual();
     1518    for (i = 0; i < 3; i++) {
     1519      rowLower[i] = -1.0e20;
     1520      colUpper[i + 2] = 0.0;
     1521    }
     1522    solution.setLogLevel(3);
     1523    solution.dual();
     1524    double rowObjective[] = { 1.0, 0.5, -10.0 };
     1525    solution.loadProblem(matrix, colLower, colUpper, objective,
     1526      rowLower, rowUpper, rowObjective);
     1527    solution.dual();
     1528    solution.loadProblem(matrix, colLower, colUpper, objective,
     1529      rowLower, rowUpper, rowObjective);
     1530    solution.primal();
     1531  }
    15351532#ifndef COIN_NO_CLP_MESSAGE
    1536      {
    1537           CoinMpsIO m;
    1538           std::string fn = dirSample + "exmip1";
    1539           if (m.readMps(fn.c_str(), "mps") == 0) {
    1540                ClpSimplex solution;
    1541                solution.loadProblem(*m.getMatrixByCol(), m.getColLower(), m.getColUpper(),
    1542                                     m.getObjCoefficients(),
    1543                                     m.getRowLower(), m.getRowUpper());
    1544                solution.dual();
    1545                // Test event handling
    1546                MyEventHandler handler;
    1547                solution.passInEventHandler(&handler);
    1548                int numberRows = solution.numberRows();
    1549                // make sure values pass has something to do
    1550                for (int i = 0; i < numberRows; i++)
    1551                     solution.setRowStatus(i, ClpSimplex::basic);
    1552                solution.objective()[1]=-2.0;
    1553                solution.primal(1);
    1554                assert (solution.secondaryStatus() == 102); // Came out at end of pass
    1555           } else {
    1556                std::cerr << "Error reading exmip1 from sample data. Skipping test." << std::endl;
    1557           }
    1558      }
    1559      // Test Message handler
    1560      {
    1561           CoinMpsIO m;
    1562           std::string fn = dirSample + "exmip1";
    1563           //fn = "Test/subGams4";
    1564           if (m.readMps(fn.c_str(), "mps") == 0) {
    1565                ClpSimplex model;
    1566                model.loadProblem(*m.getMatrixByCol(), m.getColLower(), m.getColUpper(),
    1567                                  m.getObjCoefficients(),
    1568                                  m.getRowLower(), m.getRowUpper());
    1569                // Message handler
    1570                MyMessageHandler messageHandler(&model);
    1571                std::cout << "Testing derived message handler" << std::endl;
    1572                model.passInMessageHandler(&messageHandler);
    1573                model.messagesPointer()->setDetailMessage(1, 102);
    1574                model.setFactorizationFrequency(10);
    1575                model.primal();
    1576                model.primal(0, 3);
    1577                model.setObjCoeff(3, -2.9473684210526314);
    1578                model.primal(0, 3);
    1579                // Write saved solutions
    1580                int nc = model.getNumCols();
    1581                size_t s;
    1582                std::deque<StdVectorDouble> fep = messageHandler.getFeasibleExtremePoints();
    1583                size_t numSavedSolutions = fep.size();
    1584                for ( s = 0; s < numSavedSolutions; ++s ) {
    1585                     const StdVectorDouble & solnVec = fep[s];
    1586                     for ( int c = 0; c < nc; ++c ) {
    1587                          if (fabs(solnVec[c]) > 1.0e-8)
    1588                               std::cout << "Saved Solution: " << s << " ColNum: " << c << " Value: " << solnVec[c] << std::endl;
    1589                     }
    1590                }
    1591                // Solve again without scaling
    1592                // and maximize then minimize
    1593                messageHandler.clearFeasibleExtremePoints();
    1594                model.scaling(0);
    1595                model.setOptimizationDirection(-1);
    1596                model.primal();
    1597                model.setOptimizationDirection(1);
    1598                model.primal();
    1599                fep = messageHandler.getFeasibleExtremePoints();
    1600                numSavedSolutions = fep.size();
    1601                for ( s = 0; s < numSavedSolutions; ++s ) {
    1602                     const StdVectorDouble & solnVec = fep[s];
    1603                     for ( int c = 0; c < nc; ++c ) {
    1604                          if (fabs(solnVec[c]) > 1.0e-8)
    1605                               std::cout << "Saved Solution: " << s << " ColNum: " << c << " Value: " << solnVec[c] << std::endl;
    1606                     }
    1607                }
    1608           } else {
    1609                std::cerr << "Error reading exmip1 from sample data. Skipping test." << std::endl;
    1610           }
    1611      }
    1612 #endif
    1613      // Test dual ranging
    1614      {
    1615           CoinMpsIO m;
    1616           std::string fn = dirSample + "exmip1";
    1617           if (m.readMps(fn.c_str(), "mps") == 0) {
    1618                ClpSimplex model;
    1619                model.loadProblem(*m.getMatrixByCol(), m.getColLower(), m.getColUpper(),
    1620                                  m.getObjCoefficients(),
    1621                                  m.getRowLower(), m.getRowUpper());
    1622                model.primal();
    1623                int which[13] = {0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12};
    1624                double costIncrease[13];
    1625                int sequenceIncrease[13];
    1626                double costDecrease[13];
    1627                int sequenceDecrease[13];
    1628                // ranging
    1629                model.dualRanging(13, which, costIncrease, sequenceIncrease,
    1630                                  costDecrease, sequenceDecrease);
    1631                int i;
    1632                for ( i = 0; i < 13; i++)
    1633                     printf("%d increase %g %d, decrease %g %d\n",
    1634                            i, costIncrease[i], sequenceIncrease[i],
    1635                            costDecrease[i], sequenceDecrease[i]);
    1636                assert (fabs(costDecrease[3]) < 1.0e-4);
    1637                assert (fabs(costIncrease[7] - 1.0) < 1.0e-4);
    1638                model.setOptimizationDirection(-1);
    1639                {
    1640                     int j;
    1641                     double * obj = model.objective();
    1642                     int n = model.numberColumns();
    1643                     for (j = 0; j < n; j++)
    1644                          obj[j] *= -1.0;
    1645                }
    1646                double costIncrease2[13];
    1647                int sequenceIncrease2[13];
    1648                double costDecrease2[13];
    1649                int sequenceDecrease2[13];
    1650                // ranging
    1651                model.dualRanging(13, which, costIncrease2, sequenceIncrease2,
    1652                                  costDecrease2, sequenceDecrease2);
    1653                for (i = 0; i < 13; i++) {
    1654                     assert (fabs(costIncrease[i] - costDecrease2[i]) < 1.0e-6);
    1655                     assert (fabs(costDecrease[i] - costIncrease2[i]) < 1.0e-6);
    1656                     assert (sequenceIncrease[i] == sequenceDecrease2[i]);
    1657                     assert (sequenceDecrease[i] == sequenceIncrease2[i]);
    1658                }
    1659                // Now delete all rows and see what happens
    1660                model.deleteRows(model.numberRows(), which);
    1661                model.primal();
    1662                // ranging
    1663                if (!model.dualRanging(8, which, costIncrease, sequenceIncrease,
    1664                                       costDecrease, sequenceDecrease)) {
    1665                     for (i = 0; i < 8; i++) {
    1666                          printf("%d increase %g %d, decrease %g %d\n",
    1667                                 i, costIncrease[i], sequenceIncrease[i],
    1668                                 costDecrease[i], sequenceDecrease[i]);
    1669                     }
    1670                }
    1671           } else {
    1672                std::cerr << "Error reading exmip1 from sample data. Skipping test." << std::endl;
    1673           }
    1674      }
    1675      // Test primal ranging
    1676      {
    1677           CoinMpsIO m;
    1678           std::string fn = dirSample + "exmip1";
    1679           if (m.readMps(fn.c_str(), "mps") == 0) {
    1680                ClpSimplex model;
    1681                model.loadProblem(*m.getMatrixByCol(), m.getColLower(), m.getColUpper(),
    1682                                  m.getObjCoefficients(),
    1683                                  m.getRowLower(), m.getRowUpper());
    1684                model.primal();
    1685                int which[13] = {0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12};
    1686                double valueIncrease[13];
    1687                int sequenceIncrease[13];
    1688                double valueDecrease[13];
    1689                int sequenceDecrease[13];
    1690                // ranging
    1691                model.primalRanging(13, which, valueIncrease, sequenceIncrease,
    1692                                    valueDecrease, sequenceDecrease);
    1693                int i;
    1694                for ( i = 0; i < 13; i++)
    1695                     printf("%d increase %g %d, decrease %g %d\n",
    1696                            i, valueIncrease[i], sequenceIncrease[i],
    1697                            valueDecrease[i], sequenceDecrease[i]);
    1698                assert (fabs(valueIncrease[3] - 0.642857) < 1.0e-4);
    1699                assert (fabs(valueIncrease[8] - 2.95113) < 1.0e-4);
    1700           } else {
    1701                std::cerr << "Error reading exmip1 from sample data. Skipping test." << std::endl;
    1702           }
     1533  {
     1534    CoinMpsIO m;
     1535    std::string fn = dirSample + "exmip1";
     1536    if (m.readMps(fn.c_str(), "mps") == 0) {
     1537      ClpSimplex solution;
     1538      solution.loadProblem(*m.getMatrixByCol(), m.getColLower(), m.getColUpper(),
     1539        m.getObjCoefficients(),
     1540        m.getRowLower(), m.getRowUpper());
     1541      solution.dual();
     1542      // Test event handling
     1543      MyEventHandler handler;
     1544      solution.passInEventHandler(&handler);
     1545      int numberRows = solution.numberRows();
     1546      // make sure values pass has something to do
     1547      for (int i = 0; i < numberRows; i++)
     1548        solution.setRowStatus(i, ClpSimplex::basic);
     1549      solution.objective()[1] = -2.0;
     1550      solution.primal(1);
     1551      assert(solution.secondaryStatus() == 102); // Came out at end of pass
     1552    } else {
     1553      std::cerr << "Error reading exmip1 from sample data. Skipping test." << std::endl;
     1554    }
     1555  }
     1556  // Test Message handler
     1557  {
     1558    CoinMpsIO m;
     1559    std::string fn = dirSample + "exmip1";
     1560    //fn = "Test/subGams4";
     1561    if (m.readMps(fn.c_str(), "mps") == 0) {
     1562      ClpSimplex model;
     1563      model.loadProblem(*m.getMatrixByCol(), m.getColLower(), m.getColUpper(),
     1564        m.getObjCoefficients(),
     1565        m.getRowLower(), m.getRowUpper());
     1566      // Message handler
     1567      MyMessageHandler messageHandler(&model);
     1568      std::cout << "Testing derived message handler" << std::endl;
     1569      model.passInMessageHandler(&messageHandler);
     1570      model.messagesPointer()->setDetailMessage(1, 102);
     1571      model.setFactorizationFrequency(10);
     1572      model.primal();
     1573      model.primal(0, 3);
     1574      model.setObjCoeff(3, -2.9473684210526314);
     1575      model.primal(0, 3);
     1576      // Write saved solutions
     1577      int nc = model.getNumCols();
     1578      size_t s;
     1579      std::deque< StdVectorDouble > fep = messageHandler.getFeasibleExtremePoints();
     1580      size_t numSavedSolutions = fep.size();
     1581      for (s = 0; s < numSavedSolutions; ++s) {
     1582        const StdVectorDouble &solnVec = fep[s];
     1583        for (int c = 0; c < nc; ++c) {
     1584          if (fabs(solnVec[c]) > 1.0e-8)
     1585            std::cout << "Saved Solution: " << s << " ColNum: " << c << " Value: " << solnVec[c] << std::endl;
     1586        }
     1587      }
     1588      // Solve again without scaling
     1589      // and maximize then minimize
     1590      messageHandler.clearFeasibleExtremePoints();
     1591      model.scaling(0);
     1592      model.setOptimizationDirection(-1);
     1593      model.primal();
     1594      model.setOptimizationDirection(1);
     1595      model.primal();
     1596      fep = messageHandler.getFeasibleExtremePoints();
     1597      numSavedSolutions = fep.size();
     1598      for (s = 0; s < numSavedSolutions; ++s) {
     1599        const StdVectorDouble &solnVec = fep[s];
     1600        for (int c = 0; c < nc; ++c) {
     1601          if (fabs(solnVec[c]) > 1.0e-8)
     1602            std::cout << "Saved Solution: " << s << " ColNum: " << c << " Value: " << solnVec[c] << std::endl;
     1603        }
     1604      }
     1605    } else {
     1606      std::cerr << "Error reading exmip1 from sample data. Skipping test." << std::endl;
     1607    }
     1608  }
     1609#endif
     1610  // Test dual ranging
     1611  {
     1612    CoinMpsIO m;
     1613    std::string fn = dirSample + "exmip1";
     1614    if (m.readMps(fn.c_str(), "mps") == 0) {
     1615      ClpSimplex model;
     1616      model.loadProblem(*m.getMatrixByCol(), m.getColLower(), m.getColUpper(),
     1617        m.getObjCoefficients(),
     1618        m.getRowLower(), m.getRowUpper());
     1619      model.primal();
     1620      int which[13] = { 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12 };
     1621      double costIncrease[13];
     1622      int sequenceIncrease[13];
     1623      double costDecrease[13];
     1624      int sequenceDecrease[13];
     1625      // ranging
     1626      model.dualRanging(13, which, costIncrease, sequenceIncrease,
     1627        costDecrease, sequenceDecrease);
     1628      int i;
     1629      for (i = 0; i < 13; i++)
     1630        printf("%d increase %g %d, decrease %g %d\n",
     1631          i, costIncrease[i], sequenceIncrease[i],
     1632          costDecrease[i], sequenceDecrease[i]);
     1633      assert(fabs(costDecrease[3]) < 1.0e-4);
     1634      assert(fabs(costIncrease[7] - 1.0) < 1.0e-4);
     1635      model.setOptimizationDirection(-1);
     1636      {
     1637        int j;
     1638        double *obj = model.objective();
     1639        int n = model.numberColumns();
     1640        for (j = 0; j < n; j++)
     1641          obj[j] *= -1.0;
     1642      }
     1643      double costIncrease2[13];
     1644      int sequenceIncrease2[13];
     1645      double costDecrease2[13];
     1646      int sequenceDecrease2[13];
     1647      // ranging
     1648      model.dualRanging(13, which, costIncrease2, sequenceIncrease2,
     1649        costDecrease2, sequenceDecrease2);
     1650      for (i = 0; i < 13; i++) {
     1651        assert(fabs(costIncrease[i] - costDecrease2[i]) < 1.0e-6);
     1652        assert(fabs(costDecrease[i] - costIncrease2[i]) < 1.0e-6);
     1653        assert(sequenceIncrease[i] == sequenceDecrease2[i]);
     1654        assert(sequenceDecrease[i] == sequenceIncrease2[i]);
     1655      }
     1656      // Now delete all rows and see what happens
     1657      model.deleteRows(model.numberRows(), which);
     1658      model.primal();
     1659      // ranging
     1660      if (!model.dualRanging(8, which, costIncrease, sequenceIncrease,
     1661            costDecrease, sequenceDecrease)) {
     1662        for (i = 0; i < 8; i++) {
     1663          printf("%d increase %g %d, decrease %g %d\n",
     1664            i, costIncrease[i], sequenceIncrease[i],
     1665            costDecrease[i], sequenceDecrease[i]);
     1666        }
     1667      }
     1668    } else {
     1669      std::cerr << "Error reading exmip1 from sample data. Skipping test." << std::endl;
     1670    }
     1671  }
     1672  // Test primal ranging
     1673  {
     1674    CoinMpsIO m;
     1675    std::string fn = dirSample + "exmip1";
     1676    if (m.readMps(fn.c_str(), "mps") == 0) {
     1677      ClpSimplex model;
     1678      model.loadProblem(*m.getMatrixByCol(), m.getColLower(), m.getColUpper(),
     1679        m.getObjCoefficients(),
     1680        m.getRowLower(), m.getRowUpper());
     1681      model.primal();
     1682      int which[13] = { 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12 };
     1683      double valueIncrease[13];
     1684      int sequenceIncrease[13];
     1685      double valueDecrease[13];
     1686      int sequenceDecrease[13];
     1687      // ranging
     1688      model.primalRanging(13, which, valueIncrease, sequenceIncrease,
     1689        valueDecrease, sequenceDecrease);
     1690      int i;
     1691      for (i = 0; i < 13; i++)
     1692        printf("%d increase %g %d, decrease %g %d\n",
     1693          i, valueIncrease[i], sequenceIncrease[i],
     1694          valueDecrease[i], sequenceDecrease[i]);
     1695      assert(fabs(valueIncrease[3] - 0.642857) < 1.0e-4);
     1696      assert(fabs(valueIncrease[8] - 2.95113) < 1.0e-4);
     1697    } else {
     1698      std::cerr << "Error reading exmip1 from sample data. Skipping test." << std::endl;
     1699    }
    17031700#if 0
    17041701          // out until I find optimization bug
     
    17121709                              NULL, NULL, rhs, rhs, NULL);
    17131710#endif
    1714      }
    1715      // Test binv etc
    1716      {
    1717           /*
     1711  }
     1712  // Test binv etc
     1713  {
     1714    /*
    17181715             Wolsey : Page 130
    17191716             max 4x1 -  x2
     
    17261723          */
    17271724
    1728           int n_cols = 2;
    1729           int n_rows = 3;
    1730 
    1731           double obj[2] = { -4.0, 1.0};
    1732           double collb[2] = {0.0, 0.0};
    1733           double colub[2] = {COIN_DBL_MAX, COIN_DBL_MAX};
    1734           double rowlb[3] = { -COIN_DBL_MAX, -COIN_DBL_MAX, -COIN_DBL_MAX};
    1735           double rowub[3] = {14.0, 3.0, 3.0};
    1736 
    1737           int rowIndices[5] =  {0,     2,    0,    1,    2};
    1738           int colIndices[5] =  {0,     0,    1,    1,    1};
    1739           double elements[5] = {7.0, 2.0, -2.0,  1.0, -2.0};
    1740           CoinPackedMatrix M(true, rowIndices, colIndices, elements, 5);
    1741 
    1742           ClpSimplex model;
    1743           model.loadProblem(M, collb, colub, obj, rowlb, rowub);
    1744           model.dual(0, 1); // keep factorization
    1745 
    1746           //check that the tableau matches wolsey (B-1 A)
    1747           // slacks in second part of binvA
    1748           double * binvA = reinterpret_cast<double*> (malloc((n_cols + n_rows) * sizeof(double)));
    1749 
    1750           printf("B-1 A by row\n");
    1751           int i;
    1752           for( i = 0; i < n_rows; i++) {
    1753                model.getBInvARow(i, binvA, binvA + n_cols);
    1754                printf("row: %d -> ", i);
    1755                for(int j = 0; j < n_cols + n_rows; j++) {
    1756                     printf("%g, ", binvA[j]);
    1757                }
    1758                printf("\n");
    1759           }
    1760           // See if can re-use factorization AND arrays
    1761           model.primal(0, 3 + 4); // keep factorization
    1762           // And do by column
    1763           printf("B-1 A by column\n");
    1764           for( i = 0; i < n_rows + n_cols; i++) {
    1765                model.getBInvACol(i, binvA);
    1766                printf("column: %d -> ", i);
    1767                for(int j = 0; j < n_rows; j++) {
    1768                     printf("%g, ", binvA[j]);
    1769                }
    1770                printf("\n");
    1771           }
    1772           /* Do twice -
     1725    int n_cols = 2;
     1726    int n_rows = 3;
     1727
     1728    double obj[2] = { -4.0, 1.0 };
     1729    double collb[2] = { 0.0, 0.0 };
     1730    double colub[2] = { COIN_DBL_MAX, COIN_DBL_MAX };
     1731    double rowlb[3] = { -COIN_DBL_MAX, -COIN_DBL_MAX, -COIN_DBL_MAX };
     1732    double rowub[3] = { 14.0, 3.0, 3.0 };
     1733
     1734    int rowIndices[5] = { 0, 2, 0, 1, 2 };
     1735    int colIndices[5] = { 0, 0, 1, 1, 1 };
     1736    double elements[5] = { 7.0, 2.0, -2.0, 1.0, -2.0 };
     1737    CoinPackedMatrix M(true, rowIndices, colIndices, elements, 5);
     1738
     1739    ClpSimplex model;
     1740    model.loadProblem(M, collb, colub, obj, rowlb, rowub);
     1741    model.dual(0, 1); // keep factorization
     1742
     1743    //check that the tableau matches wolsey (B-1 A)
     1744    // slacks in second part of binvA
     1745    double *binvA = reinterpret_cast< double * >(malloc((n_cols + n_rows) * sizeof(double)));
     1746
     1747    printf("B-1 A by row\n");
     1748    int i;
     1749    for (i = 0; i < n_rows; i++) {
     1750      model.getBInvARow(i, binvA, binvA + n_cols);
     1751      printf("row: %d -> ", i);
     1752      for (int j = 0; j < n_cols + n_rows; j++) {
     1753        printf("%g, ", binvA[j]);
     1754      }
     1755      printf("\n");
     1756    }
     1757    // See if can re-use factorization AND arrays
     1758    model.primal(0, 3 + 4); // keep factorization
     1759    // And do by column
     1760    printf("B-1 A by column\n");
     1761    for (i = 0; i < n_rows + n_cols; i++) {
     1762      model.getBInvACol(i, binvA);
     1763      printf("column: %d -> ", i);
     1764      for (int j = 0; j < n_rows; j++) {
     1765        printf("%g, ", binvA[j]);
     1766      }
     1767      printf("\n");
     1768    }
     1769    /* Do twice -
    17731770             without and with scaling
    17741771          */
    1775           // set scaling off
    1776           model.scaling(0);
    1777           for (int iPass = 0; iPass < 2; iPass++) {
    1778                model.primal(0, 3 + 4); // keep factorization
    1779                const double * rowScale = model.rowScale();
    1780                const double * columnScale = model.columnScale();
    1781                if (!iPass)
    1782                     assert (!rowScale);
    1783                else
    1784                     assert (rowScale); // only true for this example
    1785                /* has to be exactly correct as in OsiClpsolverInterface.cpp
     1772    // set scaling off
     1773    model.scaling(0);
     1774    for (int iPass = 0; iPass < 2; iPass++) {
     1775      model.primal(0, 3 + 4); // keep factorization
     1776      const double *rowScale = model.rowScale();
     1777      const double *columnScale = model.columnScale();
     1778      if (!iPass)
     1779        assert(!rowScale);
     1780      else
     1781        assert(rowScale); // only true for this example
     1782      /* has to be exactly correct as in OsiClpsolverInterface.cpp
    17861783                  (also redo each pass as may change
    17871784               */
    1788                printf("B-1 A");
    1789                for( i = 0; i < n_rows; i++) {
    1790                     model.getBInvARow(i, binvA, binvA + n_cols);
    1791                     printf("\nrow: %d -> ", i);
    1792                     int j;
    1793                     // First columns
    1794                     for(j = 0; j < n_cols; j++) {
    1795                          if (binvA[j]) {
    1796                               printf("(%d %g), ", j, binvA[j]);
    1797                          }
    1798                     }
    1799                     // now rows
    1800                     for(j = 0; j < n_rows; j++) {
    1801                          if (binvA[j+n_cols]) {
    1802                               printf("(%d %g), ", j + n_cols, binvA[j+n_cols]);
    1803                          }
    1804                     }
    1805                }
    1806                printf("\n");
    1807                printf("And by column (trickier)");
    1808                const int * pivotVariable = model.pivotVariable();
    1809                for( i = 0; i < n_cols + n_rows; i++) {
    1810                     model.getBInvACol(i, binvA);
    1811                     printf("\ncolumn: %d -> ", i);
    1812                     for(int j = 0; j < n_rows; j++) {
    1813                          if (binvA[j]) {
    1814                               // need to know pivot variable for +1/-1 (slack) and row/column scaling
    1815                               int pivot = pivotVariable[j];
    1816                               if (pivot < n_cols) {
    1817                                    // scaled coding is in just in case
    1818                                    if (!columnScale) {
    1819                                         printf("(%d %g), ", j, binvA[j]);
    1820                                    } else {
    1821                                         printf("(%d %g), ", j, binvA[j]*columnScale[pivot]);
    1822                                    }
    1823                               } else {
    1824                                    if (!rowScale) {
    1825                                         printf("(%d %g), ", j, binvA[j]);
    1826                                    } else {
    1827                                         printf("(%d %g), ", j, binvA[j] / rowScale[pivot-n_cols]);
    1828                                    }
    1829                               }
    1830                          }
    1831                     }
    1832                }
    1833                printf("\n");
    1834                printf("binvrow");
    1835                for( i = 0; i < n_rows; i++) {
    1836                     model.getBInvRow(i, binvA);
    1837                     printf("\nrow: %d -> ", i);
    1838                     int j;
    1839                     for (j = 0; j < n_rows; j++) {
    1840                          if (binvA[j])
    1841                               printf("(%d %g), ", j, binvA[j]);
    1842                     }
    1843                }
    1844                printf("\n");
    1845                printf("And by column ");
    1846                for( i = 0; i < n_rows; i++) {
    1847                     model.getBInvCol(i, binvA);
    1848                     printf("\ncol: %d -> ", i);
    1849                     int j;
    1850                     for (j = 0; j < n_rows; j++) {
    1851                          if (binvA[j])
    1852                               printf("(%d %g), ", j, binvA[j]);
    1853                     }
    1854                }
    1855                printf("\n");
    1856                // now deal with next pass
    1857                if (!iPass) {
    1858                     // get scaling for testing
    1859                     model.scaling(1);
    1860                }
     1785      printf("B-1 A");
     1786      for (i = 0; i < n_rows; i++) {
     1787        model.getBInvARow(i, binvA, binvA + n_cols);
     1788        printf("\nrow: %d -> ", i);
     1789        int j;
     1790        // First columns
     1791        for (j = 0; j < n_cols; j++) {
     1792          if (binvA[j]) {
     1793            printf("(%d %g), ", j, binvA[j]);
    18611794          }
    1862           free(binvA);
    1863           model.setColUpper(1, 2.0);
    1864           model.dual(0, 2 + 4); // use factorization and arrays
    1865           model.dual(0, 2); // hopefully will not use factorization
    1866           model.primal(0, 3 + 4); // keep factorization
    1867           // but say basis has changed
    1868           model.setWhatsChanged(model.whatsChanged()&(~512));
    1869           model.dual(0, 2); // hopefully will not use factorization
    1870      }
    1871      // test steepest edge
    1872      {
    1873           CoinMpsIO m;
    1874           std::string fn = dirSample + "finnis";
    1875           int returnCode = m.readMps(fn.c_str(), "mps");
    1876           if (returnCode) {
    1877                // probable cause is that gz not there
    1878                fprintf(stderr, "Unable to open finnis.mps in %s\n", dirSample.c_str());
    1879                fprintf(stderr, "Most probable cause is that sample data is not available, or finnis.mps is gzipped i.e. finnis.mps.gz and libz has not been activated\n");
    1880                fprintf(stderr, "Either gunzip files or edit Makefiles/Makefile.location to get libz\n");
     1795        }
     1796        // now rows
     1797        for (j = 0; j < n_rows; j++) {
     1798          if (binvA[j + n_cols]) {
     1799            printf("(%d %g), ", j + n_cols, binvA[j + n_cols]);
     1800          }
     1801        }
     1802      }
     1803      printf("\n");
     1804      printf("And by column (trickier)");
     1805      const int *pivotVariable = model.pivotVariable();
     1806      for (i = 0; i < n_cols + n_rows; i++) {
     1807        model.getBInvACol(i, binvA);
     1808        printf("\ncolumn: %d -> ", i);
     1809        for (int j = 0; j < n_rows; j++) {
     1810          if (binvA[j]) {
     1811            // need to know pivot variable for +1/-1 (slack) and row/column scaling
     1812            int pivot = pivotVariable[j];
     1813            if (pivot < n_cols) {
     1814              // scaled coding is in just in case
     1815              if (!columnScale) {
     1816                printf("(%d %g), ", j, binvA[j]);
     1817              } else {
     1818                printf("(%d %g), ", j, binvA[j] * columnScale[pivot]);
     1819              }
     1820            } else {
     1821              if (!rowScale) {
     1822                printf("(%d %g), ", j, binvA[j]);
     1823              } else {
     1824                printf("(%d %g), ", j, binvA[j] / rowScale[pivot - n_cols]);
     1825              }
     1826            }
     1827          }
     1828        }
     1829      }
     1830      printf("\n");
     1831      printf("binvrow");
     1832      for (i = 0; i < n_rows; i++) {
     1833        model.getBInvRow(i, binvA);
     1834        printf("\nrow: %d -> ", i);
     1835        int j;
     1836        for (j = 0; j < n_rows; j++) {
     1837          if (binvA[j])
     1838            printf("(%d %g), ", j, binvA[j]);
     1839        }
     1840      }
     1841      printf("\n");
     1842      printf("And by column ");
     1843      for (i = 0; i < n_rows; i++) {
     1844        model.getBInvCol(i, binvA);
     1845        printf("\ncol: %d -> ", i);
     1846        int j;
     1847        for (j = 0; j < n_rows; j++) {
     1848          if (binvA[j])
     1849            printf("(%d %g), ", j, binvA[j]);
     1850        }
     1851      }
     1852      printf("\n");
     1853      // now deal with next pass
     1854      if (!iPass) {
     1855        // get scaling for testing
     1856        model.scaling(1);
     1857      }
     1858    }
     1859    free(binvA);
     1860    model.setColUpper(1, 2.0);
     1861    model.dual(0, 2 + 4); // use factorization and arrays
     1862    model.dual(0, 2); // hopefully will not use factorization
     1863    model.primal(0, 3 + 4); // keep factorization
     1864    // but say basis has changed
     1865    model.setWhatsChanged(model.whatsChanged() & (~512));
     1866    model.dual(0, 2); // hopefully will not use factorization
     1867  }
     1868  // test steepest edge
     1869  {
     1870    CoinMpsIO m;
     1871    std::string fn = dirSample + "finnis";
     1872    int returnCode = m.readMps(fn.c_str(), "mps");
     1873    if (returnCode) {
     1874      // probable cause is that gz not there
     1875      fprintf(stderr, "Unable to open finnis.mps in %s\n", dirSample.c_str());
     1876      fprintf(stderr, "Most probable cause is that sample data is not available, or finnis.mps is gzipped i.e. finnis.mps.gz and libz has not been activated\n");
     1877      fprintf(stderr, "Either gunzip files or edit Makefiles/Makefile.location to get libz\n");
     1878    } else {
     1879      ClpModel model;
     1880      model.loadProblem(*m.getMatrixByCol(), m.getColLower(),
     1881        m.getColUpper(),
     1882        m.getObjCoefficients(),
     1883        m.getRowLower(), m.getRowUpper());
     1884      ClpSimplex solution(model);
     1885
     1886      solution.scaling(1);
     1887      solution.setDualBound(1.0e8);
     1888      //solution.factorization()->maximumPivots(1);
     1889      //solution.setLogLevel(3);
     1890      solution.setDualTolerance(1.0e-7);
     1891      // set objective sense,
     1892      ClpDualRowSteepest steep;
     1893      solution.setDualRowPivotAlgorithm(steep);
     1894      solution.setDblParam(ClpObjOffset, m.objectiveOffset());
     1895      solution.dual();
     1896    }
     1897  }
     1898  // test normal solution
     1899  {
     1900    CoinMpsIO m;
     1901    std::string fn = dirSample + "afiro";
     1902    if (m.readMps(fn.c_str(), "mps") == 0) {
     1903      ClpSimplex solution;
     1904      ClpModel model;
     1905      // do twice - without and with scaling
     1906      int iPass;
     1907      for (iPass = 0; iPass < 2; iPass++) {
     1908        // explicit row objective for testing
     1909        int nr = m.getNumRows();
     1910        double *rowObj = new double[nr];
     1911        CoinFillN(rowObj, nr, 0.0);
     1912        model.loadProblem(*m.getMatrixByCol(), m.getColLower(), m.getColUpper(),
     1913          m.getObjCoefficients(),
     1914          m.getRowLower(), m.getRowUpper(), rowObj);
     1915        delete[] rowObj;
     1916        solution = ClpSimplex(model);
     1917        if (iPass) {
     1918          solution.scaling();
     1919        }
     1920        solution.dual();
     1921        solution.dual();
     1922        // test optimal
     1923        assert(solution.status() == 0);
     1924        int numberColumns = solution.numberColumns();
     1925        int numberRows = solution.numberRows();
     1926        CoinPackedVector colsol(numberColumns, solution.primalColumnSolution());
     1927        double *objective = solution.objective();
     1928#ifndef NDEBUG
     1929        double objValue = colsol.dotProduct(objective);
     1930#endif
     1931        CoinRelFltEq eq(1.0e-8);
     1932        assert(eq(objValue, -4.6475314286e+02));
     1933        solution.dual();
     1934        assert(eq(solution.objectiveValue(), -4.6475314286e+02));
     1935        double *lower = solution.columnLower();
     1936        double *upper = solution.columnUpper();
     1937        double *sol = solution.primalColumnSolution();
     1938        double *result = new double[numberColumns];
     1939        CoinFillN(result, numberColumns, 0.0);
     1940        solution.matrix()->transposeTimes(solution.dualRowSolution(), result);
     1941        int iRow, iColumn;
     1942        // see if feasible and dual feasible
     1943        for (iColumn = 0; iColumn < numberColumns; iColumn++) {
     1944          double value = sol[iColumn];
     1945          assert(value < upper[iColumn] + 1.0e-8);
     1946          assert(value > lower[iColumn] - 1.0e-8);
     1947          value = objective[iColumn] - result[iColumn];
     1948          assert(value > -1.0e-5);
     1949          if (sol[iColumn] > 1.0e-5)
     1950            assert(value < 1.0e-5);
     1951        }
     1952        delete[] result;
     1953        result = new double[numberRows];
     1954        CoinFillN(result, numberRows, 0.0);
     1955        solution.matrix()->times(colsol, result);
     1956        lower = solution.rowLower();
     1957        upper = solution.rowUpper();
     1958        sol = solution.primalRowSolution();
     1959#ifndef NDEBUG
     1960        for (iRow = 0; iRow < numberRows; iRow++) {
     1961          double value = result[iRow];
     1962          assert(eq(value, sol[iRow]));
     1963          assert(value < upper[iRow] + 1.0e-8);
     1964          assert(value > lower[iRow] - 1.0e-8);
     1965        }
     1966#endif
     1967        delete[] result;
     1968        // test row objective
     1969        double *rowObjective = solution.rowObjective();
     1970        CoinDisjointCopyN(solution.dualRowSolution(), numberRows, rowObjective);
     1971        CoinDisjointCopyN(solution.dualColumnSolution(), numberColumns, objective);
     1972        // this sets up all slack basis
     1973        solution.createStatus();
     1974        solution.dual();
     1975        CoinFillN(rowObjective, numberRows, 0.0);
     1976        CoinDisjointCopyN(m.getObjCoefficients(), numberColumns, objective);
     1977        solution.dual();
     1978      }
     1979    } else {
     1980      std::cerr << "Error reading afiro from sample data. Skipping test." << std::endl;
     1981    }
     1982  }
     1983  // test unbounded
     1984  {
     1985    CoinMpsIO m;
     1986    std::string fn = dirSample + "brandy";
     1987    if (m.readMps(fn.c_str(), "mps") == 0) {
     1988      ClpSimplex solution;
     1989      // do twice - without and with scaling
     1990      int iPass;
     1991      for (iPass = 0; iPass < 2; iPass++) {
     1992        solution.loadProblem(*m.getMatrixByCol(), m.getColLower(), m.getColUpper(),
     1993          m.getObjCoefficients(),
     1994          m.getRowLower(), m.getRowUpper());
     1995        if (iPass)
     1996          solution.scaling();
     1997        solution.setOptimizationDirection(-1);
     1998        // test unbounded and ray
     1999#ifdef DUAL
     2000        solution.setDualBound(100.0);
     2001        solution.dual();
     2002#else
     2003        solution.primal();
     2004#endif
     2005        assert(solution.status() == 2);
     2006        int numberColumns = solution.numberColumns();
     2007        int numberRows = solution.numberRows();
     2008        double *lower = solution.columnLower();
     2009        double *upper = solution.columnUpper();
     2010        double *sol = solution.primalColumnSolution();
     2011        double *ray = solution.unboundedRay();
     2012        double *objective = solution.objective();
     2013        double objChange = 0.0;
     2014        int iRow, iColumn;
     2015        // make sure feasible and columns form ray
     2016        for (iColumn = 0; iColumn < numberColumns; iColumn++) {
     2017          double value = sol[iColumn];
     2018          assert(value < upper[iColumn] + 1.0e-8);
     2019          assert(value > lower[iColumn] - 1.0e-8);
     2020          value = ray[iColumn];
     2021          if (value > 0.0)
     2022            assert(upper[iColumn] > 1.0e30);
     2023          else if (value < 0.0)
     2024            assert(lower[iColumn] < -1.0e30);
     2025          objChange += value * objective[iColumn];
     2026        }
     2027        // make sure increasing objective
     2028        assert(objChange > 0.0);
     2029        double *result = new double[numberRows];
     2030        CoinFillN(result, numberRows, 0.0);
     2031        solution.matrix()->times(sol, result);
     2032        lower = solution.rowLower();
     2033        upper = solution.rowUpper();
     2034        sol = solution.primalRowSolution();
     2035#ifndef NDEBUG
     2036        for (iRow = 0; iRow < numberRows; iRow++) {
     2037          double value = result[iRow];
     2038          assert(eq(value, sol[iRow]));
     2039          assert(value < upper[iRow] + 2.0e-8);
     2040          assert(value > lower[iRow] - 2.0e-8);
     2041        }
     2042#endif
     2043        CoinFillN(result, numberRows, 0.0);
     2044        solution.matrix()->times(ray, result);
     2045        // there may be small differences (especially if scaled)
     2046        for (iRow = 0; iRow < numberRows; iRow++) {
     2047          double value = result[iRow];
     2048          if (value > 1.0e-8)
     2049            assert(upper[iRow] > 1.0e30);
     2050          else if (value < -1.0e-8)
     2051            assert(lower[iRow] < -1.0e30);
     2052        }
     2053        delete[] result;
     2054        delete[] ray;
     2055      }
     2056    } else {
     2057      std::cerr << "Error reading brandy from sample data. Skipping test." << std::endl;
     2058    }
     2059  }
     2060  // test infeasible
     2061  {
     2062    CoinMpsIO m;
     2063    std::string fn = dirSample + "brandy";
     2064    if (m.readMps(fn.c_str(), "mps") == 0) {
     2065      ClpSimplex solution;
     2066      // do twice - without and with scaling
     2067      int iPass;
     2068      for (iPass = 0; iPass < 2; iPass++) {
     2069        solution.loadProblem(*m.getMatrixByCol(), m.getColLower(), m.getColUpper(),
     2070          m.getObjCoefficients(),
     2071          m.getRowLower(), m.getRowUpper());
     2072        if (iPass)
     2073          solution.scaling();
     2074        // test infeasible and ray
     2075        solution.columnUpper()[0] = 0.0;
     2076#if 1 //def DUAL
     2077        solution.setDualBound(100.0);
     2078        solution.dual();
     2079#else
     2080        solution.primal();
     2081#endif
     2082        assert(solution.status() == 1);
     2083        int numberColumns = solution.numberColumns();
     2084        int numberRows = solution.numberRows();
     2085        double *lower = solution.rowLower();
     2086        double *upper = solution.rowUpper();
     2087        double *ray = solution.infeasibilityRay();
     2088        assert(ray);
     2089        // construct proof of infeasibility
     2090        int iRow, iColumn;
     2091        double lo = 0.0, up = 0.0;
     2092        int nl = 0, nu = 0;
     2093        for (iRow = 0; iRow < numberRows; iRow++) {
     2094          if (lower[iRow] > -1.0e20) {
     2095            lo += ray[iRow] * lower[iRow];
    18812096          } else {
    1882                ClpModel model;
    1883                model.loadProblem(*m.getMatrixByCol(), m.getColLower(),
    1884                                  m.getColUpper(),
    1885                                  m.getObjCoefficients(),
    1886                                  m.getRowLower(), m.getRowUpper());
    1887                ClpSimplex solution(model);
    1888 
    1889                solution.scaling(1);
    1890                solution.setDualBound(1.0e8);
    1891                //solution.factorization()->maximumPivots(1);
    1892                //solution.setLogLevel(3);
    1893                solution.setDualTolerance(1.0e-7);
    1894                // set objective sense,
    1895                ClpDualRowSteepest steep;
    1896                solution.setDualRowPivotAlgorithm(steep);
    1897                solution.setDblParam(ClpObjOffset, m.objectiveOffset());
    1898                solution.dual();
     2097            if (ray[iRow] > 1.0e-8)
     2098              nl++;
    18992099          }
    1900      }
    1901      // test normal solution
    1902      {
    1903           CoinMpsIO m;
    1904           std::string fn = dirSample + "afiro";
    1905           if (m.readMps(fn.c_str(), "mps") == 0) {
    1906                ClpSimplex solution;
    1907                ClpModel model;
    1908                // do twice - without and with scaling
    1909                int iPass;
    1910                for (iPass = 0; iPass < 2; iPass++) {
    1911                     // explicit row objective for testing
    1912                     int nr = m.getNumRows();
    1913                     double * rowObj = new double[nr];
    1914                     CoinFillN(rowObj, nr, 0.0);
    1915                     model.loadProblem(*m.getMatrixByCol(), m.getColLower(), m.getColUpper(),
    1916                                       m.getObjCoefficients(),
    1917                                       m.getRowLower(), m.getRowUpper(), rowObj);
    1918                     delete [] rowObj;
    1919                     solution = ClpSimplex(model);
    1920                     if (iPass) {
    1921                          solution.scaling();
    1922                     }
    1923                     solution.dual();
    1924                     solution.dual();
    1925                     // test optimal
    1926                     assert (solution.status() == 0);
    1927                     int numberColumns = solution.numberColumns();
    1928                     int numberRows = solution.numberRows();
    1929                     CoinPackedVector colsol(numberColumns, solution.primalColumnSolution());
    1930                     double * objective = solution.objective();
    1931 #ifndef NDEBUG
    1932                     double objValue = colsol.dotProduct(objective);
    1933 #endif
    1934                     CoinRelFltEq eq(1.0e-8);
    1935                     assert(eq(objValue, -4.6475314286e+02));
    1936                     solution.dual();
    1937                     assert(eq(solution.objectiveValue(), -4.6475314286e+02));
    1938                     double * lower = solution.columnLower();
    1939                     double * upper = solution.columnUpper();
    1940                     double * sol = solution.primalColumnSolution();
    1941                     double * result = new double[numberColumns];
    1942                     CoinFillN ( result, numberColumns, 0.0);
    1943                     solution.matrix()->transposeTimes(solution.dualRowSolution(), result);
    1944                     int iRow , iColumn;
    1945                     // see if feasible and dual feasible
    1946                     for (iColumn = 0; iColumn < numberColumns; iColumn++) {
    1947                          double value = sol[iColumn];
    1948                          assert(value < upper[iColumn] + 1.0e-8);
    1949                          assert(value > lower[iColumn] - 1.0e-8);
    1950                          value = objective[iColumn] - result[iColumn];
    1951                          assert (value > -1.0e-5);
    1952                          if (sol[iColumn] > 1.0e-5)
    1953                               assert (value < 1.0e-5);
    1954                     }
    1955                     delete [] result;
    1956                     result = new double[numberRows];
    1957                     CoinFillN ( result, numberRows, 0.0);
    1958                     solution.matrix()->times(colsol, result);
    1959                     lower = solution.rowLower();
    1960                     upper = solution.rowUpper();
    1961                     sol = solution.primalRowSolution();
    1962 #ifndef NDEBUG
    1963                     for (iRow = 0; iRow < numberRows; iRow++) {
    1964                          double value = result[iRow];
    1965                          assert(eq(value, sol[iRow]));
    1966                          assert(value < upper[iRow] + 1.0e-8);
    1967                          assert(value > lower[iRow] - 1.0e-8);
    1968                     }
    1969 #endif
    1970                     delete [] result;
    1971                     // test row objective
    1972                     double * rowObjective = solution.rowObjective();
    1973                     CoinDisjointCopyN(solution.dualRowSolution(), numberRows, rowObjective);
    1974                     CoinDisjointCopyN(solution.dualColumnSolution(), numberColumns, objective);
    1975                     // this sets up all slack basis
    1976                     solution.createStatus();
    1977                     solution.dual();
    1978                     CoinFillN(rowObjective, numberRows, 0.0);
    1979                     CoinDisjointCopyN(m.getObjCoefficients(), numberColumns, objective);
    1980                     solution.dual();
    1981                }
     2100          if (upper[iRow] < 1.0e20) {
     2101            up += ray[iRow] * upper[iRow];
    19822102          } else {
    1983                std::cerr << "Error reading afiro from sample data. Skipping test." << std::endl;
     2103            if (ray[iRow] > 1.0e-8)
     2104              nu++;
    19842105          }
    1985      }
    1986      // test unbounded
    1987      {
    1988           CoinMpsIO m;
    1989           std::string fn = dirSample + "brandy";
    1990           if (m.readMps(fn.c_str(), "mps") == 0) {
    1991                ClpSimplex solution;
    1992                // do twice - without and with scaling
    1993                int iPass;
    1994                for (iPass = 0; iPass < 2; iPass++) {
    1995                     solution.loadProblem(*m.getMatrixByCol(), m.getColLower(), m.getColUpper(),
    1996                                          m.getObjCoefficients(),
    1997                                          m.getRowLower(), m.getRowUpper());
    1998                     if (iPass)
    1999                          solution.scaling();
    2000                     solution.setOptimizationDirection(-1);
    2001                     // test unbounded and ray
    2002 #ifdef DUAL
    2003                     solution.setDualBound(100.0);
    2004                     solution.dual();
    2005 #else
    2006                     solution.primal();
    2007 #endif
    2008                     assert (solution.status() == 2);
    2009                     int numberColumns = solution.numberColumns();
    2010                     int numberRows = solution.numberRows();
    2011                     double * lower = solution.columnLower();
    2012                     double * upper = solution.columnUpper();
    2013                     double * sol = solution.primalColumnSolution();
    2014                     double * ray = solution.unboundedRay();
    2015                     double * objective = solution.objective();
    2016                     double objChange = 0.0;
    2017                     int iRow , iColumn;
    2018                     // make sure feasible and columns form ray
    2019                     for (iColumn = 0; iColumn < numberColumns; iColumn++) {
    2020                          double value = sol[iColumn];
    2021                          assert(value < upper[iColumn] + 1.0e-8);
    2022                          assert(value > lower[iColumn] - 1.0e-8);
    2023                          value = ray[iColumn];
    2024                          if (value > 0.0)
    2025                               assert(upper[iColumn] > 1.0e30);
    2026                          else if (value < 0.0)
    2027                               assert(lower[iColumn] < -1.0e30);
    2028                          objChange += value * objective[iColumn];
    2029                     }
    2030                     // make sure increasing objective
    2031                     assert(objChange > 0.0);
    2032                     double * result = new double[numberRows];
    2033                     CoinFillN ( result, numberRows, 0.0);
    2034                     solution.matrix()->times(sol, result);
    2035                     lower = solution.rowLower();
    2036                     upper = solution.rowUpper();
    2037                     sol = solution.primalRowSolution();
    2038 #ifndef NDEBUG
    2039                     for (iRow = 0; iRow < numberRows; iRow++) {
    2040                          double value = result[iRow];
    2041                          assert(eq(value, sol[iRow]));
    2042                          assert(value < upper[iRow] + 2.0e-8);
    2043                          assert(value > lower[iRow] - 2.0e-8);
    2044                     }
    2045 #endif
    2046                     CoinFillN ( result, numberRows, 0.0);
    2047                     solution.matrix()->times(ray, result);
    2048                     // there may be small differences (especially if scaled)
    2049                     for (iRow = 0; iRow < numberRows; iRow++) {
    2050                          double value = result[iRow];
    2051                          if (value > 1.0e-8)
    2052                               assert(upper[iRow] > 1.0e30);
    2053                          else if (value < -1.0e-8)
    2054                               assert(lower[iRow] < -1.0e30);
    2055                     }
    2056                     delete [] result;
    2057                     delete [] ray;
    2058                }
    2059           } else {
    2060                std::cerr << "Error reading brandy from sample data. Skipping test." << std::endl;
     2106        }
     2107        if (nl)
     2108          lo = -1.0e100;
     2109        if (nu)
     2110          up = 1.0e100;
     2111        double *result = new double[numberColumns];
     2112        double lo2 = 0.0, up2 = 0.0;
     2113        CoinFillN(result, numberColumns, 0.0);
     2114        solution.matrix()->transposeTimes(ray, result);
     2115        lower = solution.columnLower();
     2116        upper = solution.columnUpper();
     2117        nl = nu = 0;
     2118        for (iColumn = 0; iColumn < numberColumns; iColumn++) {
     2119          if (result[iColumn] > 1.0e-8) {
     2120            if (lower[iColumn] > -1.0e20)
     2121              lo2 += result[iColumn] * lower[iColumn];
     2122            else
     2123              nl++;
     2124            if (upper[iColumn] < 1.0e20)
     2125              up2 += result[iColumn] * upper[iColumn];
     2126            else
     2127              nu++;
     2128          } else if (result[iColumn] < -1.0e-8) {
     2129            if (lower[iColumn] > -1.0e20)
     2130              up2 += result[iColumn] * lower[iColumn];
     2131            else
     2132              nu++;
     2133            if (upper[iColumn] < 1.0e20)
     2134              lo2 += result[iColumn] * upper[iColumn];
     2135            else
     2136              nl++;
    20612137          }
    2062      }
    2063      // test infeasible
    2064      {
    2065           CoinMpsIO m;
    2066           std::string fn = dirSample + "brandy";
    2067           if (m.readMps(fn.c_str(), "mps") == 0) {
    2068                ClpSimplex solution;
    2069                // do twice - without and with scaling
    2070                int iPass;
    2071                for (iPass = 0; iPass < 2; iPass++) {
    2072                     solution.loadProblem(*m.getMatrixByCol(), m.getColLower(), m.getColUpper(),
    2073                                          m.getObjCoefficients(),
    2074                                          m.getRowLower(), m.getRowUpper());
    2075                     if (iPass)
    2076                          solution.scaling();
    2077                     // test infeasible and ray
    2078                     solution.columnUpper()[0] = 0.0;
    2079 #if 1 //def DUAL
    2080                     solution.setDualBound(100.0);
    2081                     solution.dual();
    2082 #else
    2083                     solution.primal();
    2084 #endif
    2085                     assert (solution.status() == 1);
    2086                     int numberColumns = solution.numberColumns();
    2087                     int numberRows = solution.numberRows();
    2088                     double * lower = solution.rowLower();
    2089                     double * upper = solution.rowUpper();
    2090                     double * ray = solution.infeasibilityRay();
    2091                     assert(ray);
    2092                     // construct proof of infeasibility
    2093                     int iRow , iColumn;
    2094                     double lo = 0.0, up = 0.0;
    2095                     int nl = 0, nu = 0;
    2096                     for (iRow = 0; iRow < numberRows; iRow++) {
    2097                          if (lower[iRow] > -1.0e20) {
    2098                               lo += ray[iRow] * lower[iRow];
    2099                          } else {
    2100                               if (ray[iRow] > 1.0e-8)
    2101                                    nl++;
    2102                          }
    2103                          if (upper[iRow] < 1.0e20) {
    2104                               up += ray[iRow] * upper[iRow];
    2105                          } else {
    2106                               if (ray[iRow] > 1.0e-8)
    2107                                    nu++;
    2108                          }
    2109                     }
    2110                     if (nl)
    2111                          lo = -1.0e100;
    2112                     if (nu)
    2113                          up = 1.0e100;
    2114                     double * result = new double[numberColumns];
    2115                     double lo2 = 0.0, up2 = 0.0;
    2116                     CoinFillN ( result, numberColumns, 0.0);
    2117                     solution.matrix()->transposeTimes(ray, result);
    2118                     lower = solution.columnLower();
    2119                     upper = solution.columnUpper();
    2120                     nl = nu = 0;
    2121                     for (iColumn = 0; iColumn < numberColumns; iColumn++) {
    2122                          if (result[iColumn] > 1.0e-8) {
    2123                               if (lower[iColumn] > -1.0e20)
    2124                                    lo2 += result[iColumn] * lower[iColumn];
    2125                               else
    2126                                    nl++;
    2127                               if (upper[iColumn] < 1.0e20)
    2128                                    up2 += result[iColumn] * upper[iColumn];
    2129                               else
    2130                                    nu++;
    2131                          } else if (result[iColumn] < -1.0e-8) {
    2132                               if (lower[iColumn] > -1.0e20)
    2133                                    up2 += result[iColumn] * lower[iColumn];
    2134                               else
    2135                                    nu++;
    2136                               if (upper[iColumn] < 1.0e20)
    2137                                    lo2 += result[iColumn] * upper[iColumn];
    2138                               else
    2139                                    nl++;
    2140                          }
    2141                     }
    2142                     if (nl)
    2143                          lo2 = -1.0e100;
    2144                     if (nu)
    2145                          up2 = 1.0e100;
    2146                     // make sure inconsistency
    2147                     assert(lo2 > up || up2 < lo);
    2148                     delete [] result;
    2149                     delete [] ray;
    2150                }
    2151           } else {
    2152                std::cerr << "Error reading brandy from sample data. Skipping test." << std::endl;
     2138        }
     2139        if (nl)
     2140          lo2 = -1.0e100;
     2141        if (nu)
     2142          up2 = 1.0e100;
     2143        // make sure inconsistency
     2144        assert(lo2 > up || up2 < lo);
     2145        delete[] result;
     2146        delete[] ray;
     2147      }
     2148    } else {
     2149      std::cerr << "Error reading brandy from sample data. Skipping test." << std::endl;
     2150    }
     2151  }
     2152  // test delete and add
     2153  {
     2154    CoinMpsIO m;
     2155    std::string fn = dirSample + "brandy";
     2156    if (m.readMps(fn.c_str(), "mps") == 0) {
     2157      ClpSimplex solution;
     2158      solution.loadProblem(*m.getMatrixByCol(), m.getColLower(), m.getColUpper(),
     2159        m.getObjCoefficients(),
     2160        m.getRowLower(), m.getRowUpper());
     2161      solution.dual();
     2162      CoinRelFltEq eq(1.0e-8);
     2163      assert(eq(solution.objectiveValue(), 1.5185098965e+03));
     2164
     2165      int numberColumns = solution.numberColumns();
     2166      int numberRows = solution.numberRows();
     2167      double *saveObj = new double[numberColumns];
     2168      double *saveLower = new double[numberRows + numberColumns];
     2169      double *saveUpper = new double[numberRows + numberColumns];
     2170      int *which = new int[numberRows + numberColumns];
     2171
     2172      CoinBigIndex numberElements = m.getMatrixByCol()->getNumElements();
     2173      CoinBigIndex *starts = new CoinBigIndex[numberRows + numberColumns];
     2174      int *index = new int[numberElements];
     2175      double *element = new double[numberElements];
     2176
     2177      const CoinBigIndex *startM;
     2178      const int *lengthM;
     2179      const int *indexM;
     2180      const double *elementM;
     2181
     2182      int n, nel;
     2183
     2184      // delete non basic columns
     2185      n = 0;
     2186      nel = 0;
     2187      int iRow, iColumn;
     2188      const double *lower = m.getColLower();
     2189      const double *upper = m.getColUpper();
     2190      const double *objective = m.getObjCoefficients();
     2191      startM = m.getMatrixByCol()->getVectorStarts();
     2192      lengthM = m.getMatrixByCol()->getVectorLengths();
     2193      indexM = m.getMatrixByCol()->getIndices();
     2194      elementM = m.getMatrixByCol()->getElements();
     2195      starts[0] = 0;
     2196      for (iColumn = 0; iColumn < numberColumns; iColumn++) {
     2197        if (solution.getColumnStatus(iColumn) != ClpSimplex::basic) {
     2198          saveObj[n] = objective[iColumn];
     2199          saveLower[n] = lower[iColumn];
     2200          saveUpper[n] = upper[iColumn];
     2201          CoinBigIndex j;
     2202          for (j = startM[iColumn]; j < startM[iColumn] + lengthM[iColumn]; j++) {
     2203            index[nel] = indexM[j];
     2204            element[nel++] = elementM[j];
    21532205          }
    2154      }
    2155      // test delete and add
    2156      {
    2157           CoinMpsIO m;
    2158           std::string fn = dirSample + "brandy";
    2159           if (m.readMps(fn.c_str(), "mps") == 0) {
    2160                ClpSimplex solution;
    2161                solution.loadProblem(*m.getMatrixByCol(), m.getColLower(), m.getColUpper(),
    2162                                     m.getObjCoefficients(),
    2163                                     m.getRowLower(), m.getRowUpper());
    2164                solution.dual();
    2165                CoinRelFltEq eq(1.0e-8);
    2166                assert(eq(solution.objectiveValue(), 1.5185098965e+03));
    2167 
    2168                int numberColumns = solution.numberColumns();
    2169                int numberRows = solution.numberRows();
    2170                double * saveObj = new double [numberColumns];
    2171                double * saveLower = new double[numberRows+numberColumns];
    2172                double * saveUpper = new double[numberRows+numberColumns];
    2173                int * which = new int [numberRows+numberColumns];
    2174 
    2175                CoinBigIndex numberElements = m.getMatrixByCol()->getNumElements();
    2176                CoinBigIndex * starts = new CoinBigIndex[numberRows+numberColumns];
    2177                int * index = new int[numberElements];
    2178                double * element = new double[numberElements];
    2179 
    2180                const CoinBigIndex * startM;
    2181                const int * lengthM;
    2182                const int * indexM;
    2183                const double * elementM;
    2184 
    2185                int n, nel;
    2186 
    2187                // delete non basic columns
    2188                n = 0;
    2189                nel = 0;
    2190                int iRow , iColumn;
    2191                const double * lower = m.getColLower();
    2192                const double * upper = m.getColUpper();
    2193                const double * objective = m.getObjCoefficients();
    2194                startM = m.getMatrixByCol()->getVectorStarts();
    2195                lengthM = m.getMatrixByCol()->getVectorLengths();
    2196                indexM = m.getMatrixByCol()->getIndices();
    2197                elementM = m.getMatrixByCol()->getElements();
    2198                starts[0] = 0;
    2199                for (iColumn = 0; iColumn < numberColumns; iColumn++) {
    2200                     if (solution.getColumnStatus(iColumn) != ClpSimplex::basic) {
    2201                          saveObj[n] = objective[iColumn];
    2202                          saveLower[n] = lower[iColumn];
    2203                          saveUpper[n] = upper[iColumn];
    2204                          CoinBigIndex j;
    2205                          for (j = startM[iColumn]; j < startM[iColumn] + lengthM[iColumn]; j++) {
    2206                               index[nel] = indexM[j];
    2207                               element[nel++] = elementM[j];
    2208                          }
    2209                          which[n++] = iColumn;
    2210                          starts[n] = nel;
    2211                     }
    2212                }
    2213                solution.deleteColumns(n, which);
    2214                solution.dual();
    2215                // Put back
    2216                solution.addColumns(n, saveLower, saveUpper, saveObj,
    2217                                    starts, index, element);
    2218                solution.dual();
    2219                assert(eq(solution.objectiveValue(), 1.5185098965e+03));
    2220                // Delete all columns and add back
    2221                n = 0;
    2222                nel = 0;
    2223                starts[0] = 0;
    2224                lower = m.getColLower();
    2225                upper = m.getColUpper();
    2226                objective = m.getObjCoefficients();
    2227                for (iColumn = 0; iColumn < numberColumns; iColumn++) {
    2228                     saveObj[n] = objective[iColumn];
    2229                     saveLower[n] = lower[iColumn];
    2230                     saveUpper[n] = upper[iColumn];
    2231                     CoinBigIndex j;
    2232                     for (j = startM[iColumn]; j < startM[iColumn] + lengthM[iColumn]; j++) {
    2233                          index[nel] = indexM[j];
    2234                          element[nel++] = elementM[j];
    2235                     }
    2236                     which[n++] = iColumn;
    2237                     starts[n] = nel;
    2238                }
    2239                solution.deleteColumns(n, which);
    2240                solution.dual();
    2241                // Put back
    2242                solution.addColumns(n, saveLower, saveUpper, saveObj,
    2243                                    starts, index, element);
    2244                solution.dual();
    2245                assert(eq(solution.objectiveValue(), 1.5185098965e+03));
    2246 
    2247                // reload with original
    2248                solution.loadProblem(*m.getMatrixByCol(), m.getColLower(), m.getColUpper(),
    2249                                     m.getObjCoefficients(),
    2250                                     m.getRowLower(), m.getRowUpper());
    2251                // delete half rows
    2252                n = 0;
    2253                nel = 0;
    2254                lower = m.getRowLower();
    2255                upper = m.getRowUpper();
    2256                startM = m.getMatrixByRow()->getVectorStarts();
    2257                lengthM = m.getMatrixByRow()->getVectorLengths();
    2258                indexM = m.getMatrixByRow()->getIndices();
    2259                elementM = m.getMatrixByRow()->getElements();
    2260                starts[0] = 0;
    2261                for (iRow = 0; iRow < numberRows; iRow++) {
    2262                     if ((iRow & 1) == 0) {
    2263                          saveLower[n] = lower[iRow];
    2264                          saveUpper[n] = upper[iRow];
    2265                          CoinBigIndex j;
    2266                          for (j = startM[iRow]; j < startM[iRow] + lengthM[iRow]; j++) {
    2267                               index[nel] = indexM[j];
    2268                               element[nel++] = elementM[j];
    2269                          }
    2270                          which[n++] = iRow;
    2271                          starts[n] = nel;
    2272                     }
    2273                }
    2274                solution.deleteRows(n, which);
    2275                solution.dual();
    2276                // Put back
    2277                solution.addRows(n, saveLower, saveUpper,
    2278                                 starts, index, element);
    2279                solution.dual();
    2280                assert(eq(solution.objectiveValue(), 1.5185098965e+03));
    2281                solution.writeMps("yy.mps");
    2282                // Delete all rows
    2283                n = 0;
    2284                nel = 0;
    2285                lower = m.getRowLower();
    2286                upper = m.getRowUpper();
    2287                starts[0] = 0;
    2288                for (iRow = 0; iRow < numberRows; iRow++) {
    2289                     saveLower[n] = lower[iRow];
    2290                     saveUpper[n] = upper[iRow];
    2291                     CoinBigIndex j;
    2292                     for (j = startM[iRow]; j < startM[iRow] + lengthM[iRow]; j++) {
    2293                          index[nel] = indexM[j];
    2294                          element[nel++] = elementM[j];
    2295                     }
    2296                     which[n++] = iRow;
    2297                     starts[n] = nel;
    2298                }
    2299                solution.deleteRows(n, which);
    2300                solution.dual();
    2301                // Put back
    2302                solution.addRows(n, saveLower, saveUpper,
    2303                                 starts, index, element);
    2304                solution.dual();
    2305                solution.writeMps("xx.mps");
    2306                assert(eq(solution.objectiveValue(), 1.5185098965e+03));
    2307                // Zero out status array to give some interest
    2308                memset(solution.statusArray() + numberColumns, 0, numberRows);
    2309                solution.primal(1);
    2310                assert(eq(solution.objectiveValue(), 1.5185098965e+03));
    2311                // Delete all columns and rows
    2312                n = 0;
    2313                for (iColumn = 0; iColumn < numberColumns; iColumn++) {
    2314                     which[n++] = iColumn;
    2315                     starts[n] = nel;
    2316                }
    2317                solution.deleteColumns(n, which);
    2318                n = 0;
    2319                for (iRow = 0; iRow < numberRows; iRow++) {
    2320                     which[n++] = iRow;
    2321                     starts[n] = nel;
    2322                }
    2323                solution.deleteRows(n, which);
    2324 
    2325                delete [] saveObj;
    2326                delete [] saveLower;
    2327                delete [] saveUpper;
    2328                delete [] which;
    2329                delete [] starts;
    2330                delete [] index;
    2331                delete [] element;
    2332           } else {
    2333                std::cerr << "Error reading brandy from sample data. Skipping test." << std::endl;
     2206          which[n++] = iColumn;
     2207          starts[n] = nel;
     2208        }
     2209      }
     2210      solution.deleteColumns(n, which);
     2211      solution.dual();
     2212      // Put back
     2213      solution.addColumns(n, saveLower, saveUpper, saveObj,
     2214        starts, index, element);
     2215      solution.dual();
     2216      assert(eq(solution.objectiveValue(), 1.5185098965e+03));
     2217      // Delete all columns and add back
     2218      n = 0;
     2219      nel = 0;
     2220      starts[0] = 0;
     2221      lower = m.getColLower();
     2222      upper = m.getColUpper();
     2223      objective = m.getObjCoefficients();
     2224      for (iColumn = 0; iColumn < numberColumns; iColumn++) {
     2225        saveObj[n] = objective[iColumn];
     2226        saveLower[n] = lower[iColumn];
     2227        saveUpper[n] = upper[iColumn];
     2228        CoinBigIndex j;
     2229        for (j = startM[iColumn]; j < startM[iColumn] + lengthM[iColumn]; j++) {
     2230          index[nel] = indexM[j];
     2231          element[nel++] = elementM[j];
     2232        }
     2233        which[n++] = iColumn;
     2234        starts[n] = nel;
     2235      }
     2236      solution.deleteColumns(n, which);
     2237      solution.dual();
     2238      // Put back
     2239      solution.addColumns(n, saveLower, saveUpper, saveObj,
     2240        starts, index, element);
     2241      solution.dual();
     2242      assert(eq(solution.objectiveValue(), 1.5185098965e+03));
     2243
     2244      // reload with original
     2245      solution.loadProblem(*m.getMatrixByCol(), m.getColLower(), m.getColUpper(),
     2246        m.getObjCoefficients(),
     2247        m.getRowLower(), m.getRowUpper());
     2248      // delete half rows
     2249      n = 0;
     2250      nel = 0;
     2251      lower = m.getRowLower();
     2252      upper = m.getRowUpper();
     2253      startM = m.getMatrixByRow()->getVectorStarts();
     2254      lengthM = m.getMatrixByRow()->getVectorLengths();
     2255      indexM = m.getMatrixByRow()->getIndices();
     2256      elementM = m.getMatrixByRow()->getElements();
     2257      starts[0] = 0;
     2258      for (iRow = 0; iRow < numberRows; iRow++) {
     2259        if ((iRow & 1) == 0) {
     2260          saveLower[n] = lower[iRow];
     2261          saveUpper[n] = upper[iRow];
     2262          CoinBigIndex j;
     2263          for (j = startM[iRow]; j < startM[iRow] + lengthM[iRow]; j++) {
     2264            index[nel] = indexM[j];
     2265            element[nel++] = elementM[j];
    23342266          }
    2335      }
     2267          which[n++] = iRow;
     2268          starts[n] = nel;
     2269        }
     2270      }
     2271      solution.deleteRows(n, which);
     2272      solution.dual();
     2273      // Put back
     2274      solution.addRows(n, saveLower, saveUpper,
     2275        starts, index, element);
     2276      solution.dual();
     2277      assert(eq(solution.objectiveValue(), 1.5185098965e+03));
     2278      solution.writeMps("yy.mps");
     2279      // Delete all rows
     2280      n = 0;
     2281      nel = 0;
     2282      lower = m.getRowLower();
     2283      upper = m.getRowUpper();
     2284      starts[0] = 0;
     2285      for (iRow = 0; iRow < numberRows; iRow++) {
     2286        saveLower[n] = lower[iRow];
     2287        saveUpper[n] = upper[iRow];
     2288        CoinBigIndex j;
     2289        for (j = startM[iRow]; j < startM[iRow] + lengthM[iRow]; j++) {
     2290          index[nel] = indexM[j];
     2291          element[nel++] = elementM[j];
     2292        }
     2293        which[n++] = iRow;
     2294        starts[n] = nel;
     2295      }
     2296      solution.deleteRows(n, which);
     2297      solution.dual();
     2298      // Put back
     2299      solution.addRows(n, saveLower, saveUpper,
     2300        starts, index, element);
     2301      solution.dual();
     2302      solution.writeMps("xx.mps");
     2303      assert(eq(solution.objectiveValue(), 1.5185098965e+03));
     2304      // Zero out status array to give some interest
     2305      memset(solution.statusArray() + numberColumns, 0, numberRows);
     2306      solution.primal(1);
     2307      assert(eq(solution.objectiveValue(), 1.5185098965e+03));
     2308      // Delete all columns and rows
     2309      n = 0;
     2310      for (iColumn = 0; iColumn < numberColumns; iColumn++) {
     2311        which[n++] = iColumn;
     2312        starts[n] = nel;
     2313      }
     2314      solution.deleteColumns(n, which);
     2315      n = 0;
     2316      for (iRow = 0; iRow < numberRows; iRow++) {
     2317        which[n++] = iRow;
     2318        starts[n] = nel;
     2319      }
     2320      solution.deleteRows(n, which);
     2321
     2322      delete[] saveObj;
     2323      delete[] saveLower;
     2324      delete[] saveUpper;
     2325      delete[] which;
     2326      delete[] starts;
     2327      delete[] index;
     2328      delete[] element;
     2329    } else {
     2330      std::cerr << "Error reading brandy from sample data. Skipping test." << std::endl;
     2331    }
     2332  }
    23362333#if 1
    2337      // Test barrier
    2338      {
    2339           CoinMpsIO m;
    2340           std::string fn = dirSample + "exmip1";
    2341           if (m.readMps(fn.c_str(), "mps") == 0) {
    2342                ClpInterior solution;
    2343                solution.loadProblem(*m.getMatrixByCol(), m.getColLower(), m.getColUpper(),
    2344                                     m.getObjCoefficients(),
    2345                                     m.getRowLower(), m.getRowUpper());
    2346                solution.primalDual();
    2347           } else {
    2348                std::cerr << "Error reading exmip1 from sample data. Skipping test." << std::endl;
     2334  // Test barrier
     2335  {
     2336    CoinMpsIO m;
     2337    std::string fn = dirSample + "exmip1";
     2338    if (m.readMps(fn.c_str(), "mps") == 0) {
     2339      ClpInterior solution;
     2340      solution.loadProblem(*m.getMatrixByCol(), m.getColLower(), m.getColUpper(),
     2341        m.getObjCoefficients(),
     2342        m.getRowLower(), m.getRowUpper());
     2343      solution.primalDual();
     2344    } else {
     2345      std::cerr << "Error reading exmip1 from sample data. Skipping test." << std::endl;
     2346    }
     2347  }
     2348#endif
     2349#if COIN_BIG_INDEX == 0
     2350  // test network
     2351#define QUADRATIC
     2352  if (1) {
     2353    std::string fn = dirSample + "input.130";
     2354    int numberColumns;
     2355    int numberRows;
     2356
     2357    FILE *fp = fopen(fn.c_str(), "r");
     2358    if (!fp) {
     2359      // Try in Data/Sample
     2360      fn = "Data/Sample/input.130";
     2361      fp = fopen(fn.c_str(), "r");
     2362    }
     2363    if (!fp) {
     2364      fprintf(stderr, "Unable to open file input.130 in dirSample or Data/Sample directory\n");
     2365    } else {
     2366      int problem;
     2367      char temp[100];
     2368      // read and skip
     2369      int x = fscanf(fp, "%s", temp);
     2370      if (x < 0)
     2371        throw("bad fscanf");
     2372      assert(!strcmp(temp, "BEGIN"));
     2373      x = fscanf(fp, "%*s %*s %d %d %*s %*s %d %*s", &problem, &numberRows,
     2374        &numberColumns);
     2375      if (x < 0)
     2376        throw("bad fscanf");
     2377      // scan down to SUPPLY
     2378      while (fgets(temp, 100, fp)) {
     2379        if (!strncmp(temp, "SUPPLY", 6))
     2380          break;
     2381      }
     2382      if (strncmp(temp, "SUPPLY", 6)) {
     2383        fprintf(stderr, "Unable to find SUPPLY\n");
     2384        exit(2);
     2385      }
     2386      // get space for rhs
     2387      double *lower = new double[numberRows];
     2388      double *upper = new double[numberRows];
     2389      int i;
     2390      for (i = 0; i < numberRows; i++) {
     2391        lower[i] = 0.0;
     2392        upper[i] = 0.0;
     2393      }
     2394      // ***** Remember to convert to C notation
     2395      while (fgets(temp, 100, fp)) {
     2396        int row;
     2397        int value;
     2398        if (!strncmp(temp, "ARCS", 4))
     2399          break;
     2400        sscanf(temp, "%d %d", &row, &value);
     2401        upper[row - 1] = -value;
     2402        lower[row - 1] = -value;
     2403      }
     2404      if (strncmp(temp, "ARCS", 4)) {
     2405        fprintf(stderr, "Unable to find ARCS\n");
     2406        exit(2);
     2407      }
     2408      // number of columns may be underestimate
     2409      int *head = new int[2 * numberColumns];
     2410      int *tail = new int[2 * numberColumns];
     2411      int *ub = new int[2 * numberColumns];
     2412      int *cost = new int[2 * numberColumns];
     2413      // ***** Remember to convert to C notation
     2414      numberColumns = 0;
     2415      while (fgets(temp, 100, fp)) {
     2416        int iHead;
     2417        int iTail;
     2418        int iUb;
     2419        int iCost;
     2420        if (!strncmp(temp, "DEMAND", 6))
     2421          break;
     2422        sscanf(temp, "%d %d %d %d", &iHead, &iTail, &iCost, &iUb);
     2423        iHead--;
     2424        iTail--;
     2425        head[numberColumns] = iHead;
     2426        tail[numberColumns] = iTail;
     2427        ub[numberColumns] = iUb;
     2428        cost[numberColumns] = iCost;
     2429        numberColumns++;
     2430      }
     2431      if (strncmp(temp, "DEMAND", 6)) {
     2432        fprintf(stderr, "Unable to find DEMAND\n");
     2433        exit(2);
     2434      }
     2435      // ***** Remember to convert to C notation
     2436      while (fgets(temp, 100, fp)) {
     2437        int row;
     2438        int value;
     2439        if (!strncmp(temp, "END", 3))
     2440          break;
     2441        sscanf(temp, "%d %d", &row, &value);
     2442        upper[row - 1] = value;
     2443        lower[row - 1] = value;
     2444      }
     2445      if (strncmp(temp, "END", 3)) {
     2446        fprintf(stderr, "Unable to find END\n");
     2447        exit(2);
     2448      }
     2449      printf("Problem %d has %d rows and %d columns\n", problem,
     2450        numberRows, numberColumns);
     2451      fclose(fp);
     2452      ClpSimplex model;
     2453      // now build model
     2454
     2455      double *objective = new double[numberColumns];
     2456      double *lowerColumn = new double[numberColumns];
     2457      double *upperColumn = new double[numberColumns];
     2458
     2459      double *element = new double[2 * numberColumns];
     2460      CoinBigIndex *start = new CoinBigIndex[numberColumns + 1];
     2461      int *row = new int[2 * numberColumns];
     2462      start[numberColumns] = 2 * numberColumns;
     2463      for (i = 0; i < numberColumns; i++) {
     2464        start[i] = 2 * i;
     2465        element[2 * i] = -1.0;
     2466        element[2 * i + 1] = 1.0;
     2467        row[2 * i] = head[i];
     2468        row[2 * i + 1] = tail[i];
     2469        lowerColumn[i] = 0.0;
     2470        upperColumn[i] = ub[i];
     2471        objective[i] = cost[i];
     2472      }
     2473      // Create Packed Matrix
     2474      CoinPackedMatrix matrix;
     2475      int *lengths = NULL;
     2476      matrix.assignMatrix(true, numberRows, numberColumns,
     2477        2 * numberColumns, element, row, start, lengths);
     2478      // load model
     2479      model.loadProblem(matrix,
     2480        lowerColumn, upperColumn, objective,
     2481        lower, upper);
     2482      model.factorization()->maximumPivots(200 + model.numberRows() / 100);
     2483      model.createStatus();
     2484      double time1 = CoinCpuTime();
     2485      model.dual();
     2486      std::cout << "Network problem, ClpPackedMatrix took " << CoinCpuTime() - time1 << " seconds" << std::endl;
     2487      ClpPlusMinusOneMatrix *plusMinus = new ClpPlusMinusOneMatrix(matrix);
     2488      assert(plusMinus->getIndices()); // would be zero if not +- one
     2489      //ClpPlusMinusOneMatrix *plusminus_matrix;
     2490
     2491      //plusminus_matrix = new ClpPlusMinusOneMatrix;
     2492
     2493      //plusminus_matrix->passInCopy(numberRows, numberColumns, true, plusMinus->getMutableIndices(),
     2494      //                         plusMinus->startPositive(),plusMinus->startNegative());
     2495      model.loadProblem(*plusMinus,
     2496        lowerColumn, upperColumn, objective,
     2497        lower, upper);
     2498      //model.replaceMatrix( plusminus_matrix , true);
     2499      delete plusMinus;
     2500      //model.createStatus();
     2501      //model.initialSolve();
     2502      //model.writeMps("xx.mps");
     2503
     2504      model.factorization()->maximumPivots(200 + model.numberRows() / 100);
     2505      model.createStatus();
     2506      time1 = CoinCpuTime();
     2507      model.dual();
     2508      std::cout << "Network problem, ClpPlusMinusOneMatrix took " << CoinCpuTime() - time1 << " seconds" << std::endl;
     2509      ClpNetworkMatrix network(numberColumns, head, tail);
     2510      model.loadProblem(network,
     2511        lowerColumn, upperColumn, objective,
     2512        lower, upper);
     2513
     2514      model.factorization()->maximumPivots(200 + model.numberRows() / 100);
     2515      model.createStatus();
     2516      time1 = CoinCpuTime();
     2517      model.dual();
     2518      std::cout << "Network problem, ClpNetworkMatrix took " << CoinCpuTime() - time1 << " seconds" << std::endl;
     2519      delete[] lower;
     2520      delete[] upper;
     2521      delete[] head;
     2522      delete[] tail;
     2523      delete[] ub;
     2524      delete[] cost;
     2525      delete[] objective;
     2526      delete[] lowerColumn;
     2527      delete[] upperColumn;
     2528    }
     2529  }
     2530#ifdef QUADRATIC
     2531  // Test quadratic to solve linear
     2532  if (1) {
     2533    CoinMpsIO m;
     2534    std::string fn = dirSample + "exmip1";
     2535    if (m.readMps(fn.c_str(), "mps") == 0) {
     2536      ClpSimplex solution;
     2537      solution.loadProblem(*m.getMatrixByCol(), m.getColLower(), m.getColUpper(),
     2538        m.getObjCoefficients(),
     2539        m.getRowLower(), m.getRowUpper());
     2540      //solution.dual();
     2541      // get quadratic part
     2542      int numberColumns = solution.numberColumns();
     2543      int *start = new int[numberColumns + 1];
     2544      int *column = new int[numberColumns];
     2545      double *element = new double[numberColumns];
     2546      int i;
     2547      start[0] = 0;
     2548      int n = 0;
     2549      int kk = numberColumns - 1;
     2550      int kk2 = numberColumns - 1;
     2551      for (i = 0; i < numberColumns; i++) {
     2552        if (i >= kk) {
     2553          column[n] = i;
     2554          if (i >= kk2)
     2555            element[n] = 1.0e-1;
     2556          else
     2557            element[n] = 0.0;
     2558          n++;
     2559        }
     2560        start[i + 1] = n;
     2561      }
     2562      // Load up objective
     2563      solution.loadQuadraticObjective(numberColumns, start, column, element);
     2564      delete[] start;
     2565      delete[] column;
     2566      delete[] element;
     2567      //solution.quadraticSLP(50,1.0e-4);
     2568      CoinRelFltEq eq(1.0e-4);
     2569      //assert(eq(objValue,820.0));
     2570      //solution.setLogLevel(63);
     2571      solution.primal();
     2572      printSol(solution);
     2573      //assert(eq(objValue,3.2368421));
     2574      //exit(77);
     2575    } else {
     2576      std::cerr << "Error reading exmip1 from sample data. Skipping test." << std::endl;
     2577    }
     2578  }
     2579  // Test quadratic
     2580  if (1) {
     2581    CoinMpsIO m;
     2582    std::string fn = dirSample + "share2qp";
     2583    //fn = "share2qpb";
     2584    if (m.readMps(fn.c_str(), "mps") == 0) {
     2585      ClpSimplex model;
     2586      model.loadProblem(*m.getMatrixByCol(), m.getColLower(), m.getColUpper(),
     2587        m.getObjCoefficients(),
     2588        m.getRowLower(), m.getRowUpper());
     2589      model.dual();
     2590      // get quadratic part
     2591      int *start = NULL;
     2592      int *column = NULL;
     2593      double *element = NULL;
     2594      m.readQuadraticMps(NULL, start, column, element, 2);
     2595      int column2[200];
     2596      double element2[200];
     2597      int start2[80];
     2598      int j;
     2599      start2[0] = 0;
     2600      int nel = 0;
     2601      bool good = false;
     2602      for (j = 0; j < 79; j++) {
     2603        if (start[j] == start[j + 1]) {
     2604          column2[nel] = j;
     2605          element2[nel] = 0.0;
     2606          nel++;
     2607        } else {
     2608          int i;
     2609          for (i = start[j]; i < start[j + 1]; i++) {
     2610            column2[nel] = column[i];
     2611            element2[nel++] = element[i];
    23492612          }
    2350      }
    2351 #endif
    2352 #if COIN_BIG_INDEX==0
    2353      // test network
    2354 #define QUADRATIC
    2355      if (1) {
    2356           std::string fn = dirSample + "input.130";
    2357           int numberColumns;
    2358           int numberRows;
    2359 
    2360           FILE * fp = fopen(fn.c_str(), "r");
    2361           if (!fp) {
    2362                // Try in Data/Sample
    2363                fn = "Data/Sample/input.130";
    2364                fp = fopen(fn.c_str(), "r");
    2365           }
    2366           if (!fp) {
    2367                fprintf(stderr, "Unable to open file input.130 in dirSample or Data/Sample directory\n");
    2368           } else {
    2369                int problem;
    2370                char temp[100];
    2371                // read and skip
    2372                int x = fscanf(fp, "%s", temp);
    2373                if (x < 0)
    2374                     throw("bad fscanf");
    2375                assert (!strcmp(temp, "BEGIN"));
    2376                x = fscanf(fp, "%*s %*s %d %d %*s %*s %d %*s", &problem, &numberRows,
    2377                           &numberColumns);
    2378                if (x < 0)
    2379                     throw("bad fscanf");
    2380                // scan down to SUPPLY
    2381                while (fgets(temp, 100, fp)) {
    2382                     if (!strncmp(temp, "SUPPLY", 6))
    2383                          break;
    2384                }
    2385                if (strncmp(temp, "SUPPLY", 6)) {
    2386                     fprintf(stderr, "Unable to find SUPPLY\n");
    2387                     exit(2);
    2388                }
    2389                // get space for rhs
    2390                double * lower = new double[numberRows];
    2391                double * upper = new double[numberRows];
    2392                int i;
    2393                for (i = 0; i < numberRows; i++) {
    2394                     lower[i] = 0.0;
    2395                     upper[i] = 0.0;
    2396                }
    2397                // ***** Remember to convert to C notation
    2398                while (fgets(temp, 100, fp)) {
    2399                     int row;
    2400                     int value;
    2401                     if (!strncmp(temp, "ARCS", 4))
    2402                          break;
    2403                     sscanf(temp, "%d %d", &row, &value);
    2404                     upper[row-1] = -value;
    2405                     lower[row-1] = -value;
    2406                }
    2407                if (strncmp(temp, "ARCS", 4)) {
    2408                     fprintf(stderr, "Unable to find ARCS\n");
    2409                     exit(2);
    2410                }
    2411                // number of columns may be underestimate
    2412                int * head = new int[2*numberColumns];
    2413                int * tail = new int[2*numberColumns];
    2414                int * ub = new int[2*numberColumns];
    2415                int * cost = new int[2*numberColumns];
    2416                // ***** Remember to convert to C notation
    2417                numberColumns = 0;
    2418                while (fgets(temp, 100, fp)) {
    2419                     int iHead;
    2420                     int iTail;
    2421                     int iUb;
    2422                     int iCost;
    2423                     if (!strncmp(temp, "DEMAND", 6))
    2424                          break;
    2425                     sscanf(temp, "%d %d %d %d", &iHead, &iTail, &iCost, &iUb);
    2426                     iHead--;
    2427                     iTail--;
    2428                     head[numberColumns] = iHead;
    2429                     tail[numberColumns] = iTail;
    2430                     ub[numberColumns] = iUb;
    2431                     cost[numberColumns] = iCost;
    2432                     numberColumns++;
    2433                }
    2434                if (strncmp(temp, "DEMAND", 6)) {
    2435                     fprintf(stderr, "Unable to find DEMAND\n");
    2436                     exit(2);
    2437                }
    2438                // ***** Remember to convert to C notation
    2439                while (fgets(temp, 100, fp)) {
    2440                     int row;
    2441                     int value;
    2442                     if (!strncmp(temp, "END", 3))
    2443                          break;
    2444                     sscanf(temp, "%d %d", &row, &value);
    2445                     upper[row-1] = value;
    2446                     lower[row-1] = value;
    2447                }
    2448                if (strncmp(temp, "END", 3)) {
    2449                     fprintf(stderr, "Unable to find END\n");
    2450                     exit(2);
    2451                }
    2452                printf("Problem %d has %d rows and %d columns\n", problem,
    2453                       numberRows, numberColumns);
    2454                fclose(fp);
    2455                ClpSimplex  model;
    2456                // now build model
    2457 
    2458                double * objective = new double[numberColumns];
    2459                double * lowerColumn = new double[numberColumns];
    2460                double * upperColumn = new double[numberColumns];
    2461 
    2462                double * element = new double [2*numberColumns];
    2463                CoinBigIndex * start = new CoinBigIndex [numberColumns+1];
    2464                int * row = new int[2*numberColumns];
    2465                start[numberColumns] = 2 * numberColumns;
    2466                for (i = 0; i < numberColumns; i++) {
    2467                     start[i] = 2 * i;
    2468                     element[2*i] = -1.0;
    2469                     element[2*i+1] = 1.0;
    2470                     row[2*i] = head[i];
    2471                     row[2*i+1] = tail[i];
    2472                     lowerColumn[i] = 0.0;
    2473                     upperColumn[i] = ub[i];
    2474                     objective[i] = cost[i];
    2475                }
    2476                // Create Packed Matrix
    2477                CoinPackedMatrix matrix;
    2478                int * lengths = NULL;
    2479                matrix.assignMatrix(true, numberRows, numberColumns,
    2480                                    2 * numberColumns, element, row, start, lengths);
    2481                // load model
    2482                model.loadProblem(matrix,
    2483                                  lowerColumn, upperColumn, objective,
    2484                                  lower, upper);
    2485                model.factorization()->maximumPivots(200 + model.numberRows() / 100);
    2486                model.createStatus();
    2487                double time1 = CoinCpuTime();
    2488                model.dual();
    2489                std::cout << "Network problem, ClpPackedMatrix took " << CoinCpuTime() - time1 << " seconds" << std::endl;
    2490                ClpPlusMinusOneMatrix * plusMinus = new ClpPlusMinusOneMatrix(matrix);
    2491                assert (plusMinus->getIndices()); // would be zero if not +- one
    2492                //ClpPlusMinusOneMatrix *plusminus_matrix;
    2493 
    2494                //plusminus_matrix = new ClpPlusMinusOneMatrix;
    2495 
    2496                //plusminus_matrix->passInCopy(numberRows, numberColumns, true, plusMinus->getMutableIndices(),
    2497                //                         plusMinus->startPositive(),plusMinus->startNegative());
    2498                model.loadProblem(*plusMinus,
    2499                                  lowerColumn, upperColumn, objective,
    2500                                  lower, upper);
    2501                //model.replaceMatrix( plusminus_matrix , true);
    2502                delete plusMinus;
    2503                //model.createStatus();
    2504                //model.initialSolve();
    2505                //model.writeMps("xx.mps");
    2506 
    2507                model.factorization()->maximumPivots(200 + model.numberRows() / 100);
    2508                model.createStatus();
    2509                time1 = CoinCpuTime();
    2510                model.dual();
    2511                std::cout << "Network problem, ClpPlusMinusOneMatrix took " << CoinCpuTime() - time1 << " seconds" << std::endl;
    2512                ClpNetworkMatrix network(numberColumns, head, tail);
    2513                model.loadProblem(network,
    2514                                  lowerColumn, upperColumn, objective,
    2515                                  lower, upper);
    2516 
    2517                model.factorization()->maximumPivots(200 + model.numberRows() / 100);
    2518                model.createStatus();
    2519                time1 = CoinCpuTime();
    2520                model.dual();
    2521                std::cout << "Network problem, ClpNetworkMatrix took " << CoinCpuTime() - time1 << " seconds" << std::endl;
    2522                delete [] lower;
    2523                delete [] upper;
    2524                delete [] head;
    2525                delete [] tail;
    2526                delete [] ub;
    2527                delete [] cost;
    2528                delete [] objective;
    2529                delete [] lowerColumn;
    2530                delete [] upperColumn;
    2531           }
    2532      }
    2533 #ifdef QUADRATIC
    2534      // Test quadratic to solve linear
    2535      if (1) {
    2536           CoinMpsIO m;
    2537           std::string fn = dirSample + "exmip1";
    2538           if (m.readMps(fn.c_str(), "mps") == 0) {
    2539                ClpSimplex solution;
    2540                solution.loadProblem(*m.getMatrixByCol(), m.getColLower(), m.getColUpper(),
    2541                                     m.getObjCoefficients(),
    2542                                     m.getRowLower(), m.getRowUpper());
    2543                //solution.dual();
    2544                // get quadratic part
    2545                int numberColumns = solution.numberColumns();
    2546                int * start = new int [numberColumns+1];
    2547                int * column = new int[numberColumns];
    2548                double * element = new double[numberColumns];
    2549                int i;
    2550                start[0] = 0;
    2551                int n = 0;
    2552                int kk = numberColumns - 1;
    2553                int kk2 = numberColumns - 1;
    2554                for (i = 0; i < numberColumns; i++) {
    2555                     if (i >= kk) {
    2556                          column[n] = i;
    2557                          if (i >= kk2)
    2558                               element[n] = 1.0e-1;
    2559                          else
    2560                               element[n] = 0.0;
    2561                          n++;
    2562                     }
    2563                     start[i+1] = n;
    2564                }
    2565                // Load up objective
    2566                solution.loadQuadraticObjective(numberColumns, start, column, element);
    2567                delete [] start;
    2568                delete [] column;
    2569                delete [] element;
    2570                //solution.quadraticSLP(50,1.0e-4);
    2571                CoinRelFltEq eq(1.0e-4);
    2572                //assert(eq(objValue,820.0));
    2573                //solution.setLogLevel(63);
    2574                solution.primal();
    2575                printSol(solution);
    2576                //assert(eq(objValue,3.2368421));
    2577                //exit(77);
    2578           } else {
    2579                std::cerr << "Error reading exmip1 from sample data. Skipping test." << std::endl;
    2580           }
    2581      }
    2582      // Test quadratic
    2583      if (1) {
    2584           CoinMpsIO m;
    2585           std::string fn = dirSample + "share2qp";
    2586           //fn = "share2qpb";
    2587           if (m.readMps(fn.c_str(), "mps") == 0) {
    2588                ClpSimplex model;
    2589                model.loadProblem(*m.getMatrixByCol(), m.getColLower(), m.getColUpper(),
    2590                                  m.getObjCoefficients(),
    2591                                  m.getRowLower(), m.getRowUpper());
    2592                model.dual();
    2593                // get quadratic part
    2594                int * start = NULL;
    2595                int * column = NULL;
    2596                double * element = NULL;
    2597                m.readQuadraticMps(NULL, start, column, element, 2);
    2598                int column2[200];
    2599                double element2[200];
    2600                int start2[80];
    2601                int j;
    2602                start2[0] = 0;
    2603                int nel = 0;
    2604                bool good = false;
    2605                for (j = 0; j < 79; j++) {
    2606                     if (start[j] == start[j+1]) {
    2607                          column2[nel] = j;
    2608                          element2[nel] = 0.0;
    2609                          nel++;
    2610                     } else {
    2611                          int i;
    2612                          for (i = start[j]; i < start[j+1]; i++) {
    2613                               column2[nel] = column[i];
    2614                               element2[nel++] = element[i];
    2615                          }
    2616                     }
    2617                     start2[j+1] = nel;
    2618                }
    2619                // Load up objective
    2620                if (good)
    2621                     model.loadQuadraticObjective(model.numberColumns(), start2, column2, element2);
    2622                else
    2623                     model.loadQuadraticObjective(model.numberColumns(), start, column, element);
    2624                delete [] start;
    2625                delete [] column;
    2626                delete [] element;
    2627                int numberColumns = model.numberColumns();
    2628                model.scaling(0);
     2613        }
     2614        start2[j + 1] = nel;
     2615      }
     2616      // Load up objective
     2617      if (good)
     2618        model.loadQuadraticObjective(model.numberColumns(), start2, column2, element2);
     2619      else
     2620        model.loadQuadraticObjective(model.numberColumns(), start, column, element);
     2621      delete[] start;
     2622      delete[] column;
     2623      delete[] element;
     2624      int numberColumns = model.numberColumns();
     2625      model.scaling(0);
    26292626#if 0
    26302627               model.nonlinearSLP(50, 1.0e-4);
    26312628#else
    2632                // Get feasible
    2633                ClpObjective * saveObjective = model.objectiveAsObject()->clone();
    2634                ClpLinearObjective zeroObjective(NULL, numberColumns);
    2635                model.setObjective(&zeroObjective);
    2636                model.dual();
    2637                model.setObjective(saveObjective);
    2638                delete saveObjective;
    2639 #endif
    2640                //model.setLogLevel(63);
    2641                //exit(77);
    2642                model.setFactorizationFrequency(10);
    2643                model.primal();
    2644                printSol(model);
     2629      // Get feasible
     2630      ClpObjective *saveObjective = model.objectiveAsObject()->clone();
     2631      ClpLinearObjective zeroObjective(NULL, numberColumns);
     2632      model.setObjective(&zeroObjective);
     2633      model.dual();
     2634      model.setObjective(saveObjective);
     2635      delete saveObjective;
     2636#endif
     2637      //model.setLogLevel(63);
     2638      //exit(77);
     2639      model.setFactorizationFrequency(10);
     2640      model.primal();
     2641      printSol(model);
    26452642#ifndef NDEBUG
    2646                double objValue = model.getObjValue();
    2647 #endif
    2648                CoinRelFltEq eq(1.0e-4);
    2649                assert(eq(objValue, -400.92));
    2650                // and again for barrier
    2651                model.barrier(false);
    2652                //printSol(model);
    2653                model.allSlackBasis();
    2654                model.primal();
    2655                //printSol(model);
    2656           } else {
    2657                std::cerr << "Error reading share2qp from sample data. Skipping test." << std::endl;
    2658           }
    2659      }
    2660      if (0) {
    2661           CoinMpsIO m;
    2662           std::string fn = "./beale";
    2663           //fn = "./jensen";
    2664           if (m.readMps(fn.c_str(), "mps") == 0) {
    2665                ClpSimplex solution;
    2666                solution.loadProblem(*m.getMatrixByCol(), m.getColLower(), m.getColUpper(),
    2667                                     m.getObjCoefficients(),
    2668                                     m.getRowLower(), m.getRowUpper());
    2669                solution.setDblParam(ClpObjOffset, m.objectiveOffset());
    2670                solution.dual();
    2671                // get quadratic part
    2672                int * start = NULL;
    2673                int * column = NULL;
    2674                double * element = NULL;
    2675                m.readQuadraticMps(NULL, start, column, element, 2);
    2676                // Load up objective
    2677                solution.loadQuadraticObjective(solution.numberColumns(), start, column, element);
    2678                delete [] start;
    2679                delete [] column;
    2680                delete [] element;
    2681                solution.primal(1);
    2682                solution.nonlinearSLP(50, 1.0e-4);
    2683                double objValue = solution.getObjValue();
    2684                CoinRelFltEq eq(1.0e-4);
    2685                assert(eq(objValue, 0.5));
    2686                solution.primal();
    2687                objValue = solution.getObjValue();
    2688                assert(eq(objValue, 0.5));
    2689           } else {
    2690                std::cerr << "Error reading beale.mps. Skipping test." << std::endl;
    2691           }
    2692      }
    2693 #endif
    2694      // Test CoinStructuredModel
    2695      {
    2696 
    2697           // Sub block
    2698           CoinModel sub;
    2699           {
    2700                // matrix data
    2701                //deliberate hiccup of 2 between 0 and 1
    2702                CoinBigIndex start[5] = {0, 4, 7, 8, 9};
    2703                int length[5] = {2, 3, 1, 1, 1};
    2704                int rows[11] = {0, 2, -1, -1, 0, 1, 2, 0, 1, 2};
    2705                double elements[11] = {7.0, 2.0, 1.0e10, 1.0e10, -2.0, 1.0, -2.0, 1, 1, 1};
    2706                CoinPackedMatrix matrix(true, 3, 5, 8, elements, rows, start, length);
    2707                // by row
    2708                matrix.reverseOrdering();
    2709                const double * element = matrix.getElements();
    2710                const int * column = matrix.getIndices();
    2711                const CoinBigIndex * rowStart = matrix.getVectorStarts();
    2712                const int * rowLength = matrix.getVectorLengths();
    2713 
    2714                // rim data
    2715                //double objective[5]={-4.0,1.0,0.0,0.0,0.0};
    2716                double rowLower[3] = {14.0, 3.0, 3.0};
    2717                double rowUpper[3] = {14.0, 3.0, 3.0};
    2718                //double columnLower[5]={0.0,0.0,0.0,0.0,0.0};
    2719                //double columnUpper[5]={100.0,100.0,100.0,100.0,100.0};
    2720 
    2721                for (int i = 0; i < 3; i++) {
    2722                     sub.addRow(rowLength[i], column + rowStart[i],
    2723                                element + rowStart[i], rowLower[i], rowUpper[i]);
    2724                }
    2725                //for (int i=0;i<5;i++) {
    2726                //sub.setColumnBounds(i,columnLower[i],columnUpper[i]);
    2727                //sub.setColumnObjective(i,objective[i]);
    2728                //}
    2729                sub.convertMatrix();
    2730           }
    2731           // Top
    2732           CoinModel top;
    2733           {
    2734                // matrix data
    2735                CoinBigIndex start[5] = {0, 2, 4, 6, 8};
    2736                int length[5] = {2, 2, 2, 2, 2};
    2737                int rows[10] = {0, 1, 0, 1, 0, 1, 0, 1, 0, 1};
    2738                double elements[10] = {1, 1, 1, 1, 1, 1, 1, 1, 1, 1};
    2739                CoinPackedMatrix matrix(true, 2, 5, 8, elements, rows, start, length);
    2740                // by row
    2741                matrix.reverseOrdering();
    2742                const double * element = matrix.getElements();
    2743                const int * column = matrix.getIndices();
    2744                const CoinBigIndex * rowStart = matrix.getVectorStarts();
    2745                const int * rowLength = matrix.getVectorLengths();
    2746 
    2747                // rim data
    2748                double objective[5] = { -4.0, 1.0, 0.0, 0.0, 0.0};
    2749                //double rowLower[3]={14.0,3.0,3.0};
    2750                //double rowUpper[3]={14.0,3.0,3.0};
    2751                double columnLower[5] = {0.0, 0.0, 0.0, 0.0, 0.0};
    2752                double columnUpper[5] = {100.0, 100.0, 100.0, 100.0, 100.0};
    2753 
    2754                for (int i = 0; i < 2; i++) {
    2755                     top.addRow(rowLength[i], column + rowStart[i],
    2756                                element + rowStart[i],
    2757                                -COIN_DBL_MAX, COIN_DBL_MAX);
    2758                }
    2759                for (int i = 0; i < 5; i++) {
    2760                     top.setColumnBounds(i, columnLower[i], columnUpper[i]);
    2761                     top.setColumnObjective(i, objective[i]);
    2762                }
    2763                top.convertMatrix();
    2764           }
    2765           // Create a structured model
    2766           CoinStructuredModel structured;
    2767           int numberBlocks = 5;
    2768           for (int i = 0; i < numberBlocks; i++) {
    2769                std::string topName = "row_master";
    2770                std::string blockName = "block_";
    2771                char bName = static_cast<char>('a' + static_cast<char>(i));
    2772                blockName.append(1, bName);
    2773                structured.addBlock(topName, blockName, top);
    2774                structured.addBlock(blockName, blockName, sub);
    2775           }
    2776           // Set rhs on first block
    2777           CoinModel * first = structured.coinBlock(0);
    2778           for (int i = 0; i < 2; i++) {
    2779                first->setRowLower(i, 0.0);
    2780                first->setRowUpper(i, 100.0);
    2781           }
    2782           // Refresh whats set
    2783           structured.refresh(0);
    2784           // Could perturb stuff, but for first go don't bother
    2785           ClpSimplex fullModel;
    2786           // There is no original stuff set - think
    2787           fullModel.loadProblem(structured, false);
    2788           fullModel.dual();
    2789           fullModel.dropNames();
    2790           fullModel.writeMps("test.mps");
    2791           // Make up very simple nested model - not realistic
    2792           // Create a structured model
    2793           CoinStructuredModel structured2;
    2794           numberBlocks = 3;
    2795           for (int i = 0; i < numberBlocks; i++) {
    2796                std::string blockName = "block_";
    2797                char bName = static_cast<char>('a' + static_cast<char>(i));
    2798                blockName.append(1, bName);
    2799                structured2.addBlock(blockName, blockName, structured);
    2800           }
    2801           fullModel.loadProblem(structured2, false);
    2802           fullModel.dual();
    2803           fullModel.dropNames();
    2804           fullModel.writeMps("test2.mps");
    2805      }
     2643      double objValue = model.getObjValue();
     2644#endif
     2645      CoinRelFltEq eq(1.0e-4);
     2646      assert(eq(objValue, -400.92));
     2647      // and again for barrier
     2648      model.barrier(false);
     2649      //printSol(model);
     2650      model.allSlackBasis();
     2651      model.primal();
     2652      //printSol(model);
     2653    } else {
     2654      std::cerr << "Error reading share2qp from sample data. Skipping test." << std::endl;
     2655    }
     2656  }
     2657  if (0) {
     2658    CoinMpsIO m;
     2659    std::string fn = "./beale";
     2660    //fn = "./jensen";
     2661    if (m.readMps(fn.c_str(), "mps") == 0) {
     2662      ClpSimplex solution;
     2663      solution.loadProblem(*m.getMatrixByCol(), m.getColLower(), m.getColUpper(),
     2664        m.getObjCoefficients(),
     2665        m.getRowLower(), m.getRowUpper());
     2666      solution.setDblParam(ClpObjOffset, m.objectiveOffset());
     2667      solution.dual();
     2668      // get quadratic part
     2669      int *start = NULL;
     2670      int *column = NULL;
     2671      double *element = NULL;
     2672      m.readQuadraticMps(NULL, start, column, element, 2);
     2673      // Load up objective
     2674      solution.loadQuadraticObjective(solution.numberColumns(), start, column, element);
     2675      delete[] start;
     2676      delete[] column;
     2677      delete[] element;
     2678      solution.primal(1);
     2679      solution.nonlinearSLP(50, 1.0e-4);
     2680      double objValue = solution.getObjValue();
     2681      CoinRelFltEq eq(1.0e-4);
     2682      assert(eq(objValue, 0.5));
     2683      solution.primal();
     2684      objValue = solution.getObjValue();
     2685      assert(eq(objValue, 0.5));
     2686    } else {
     2687      std::cerr << "Error reading beale.mps. Skipping test." << std::endl;
     2688    }
     2689  }
     2690#endif
     2691  // Test CoinStructuredModel
     2692  {
     2693
     2694    // Sub block
     2695    CoinModel sub;
     2696    {
     2697      // matrix data
     2698      //deliberate hiccup of 2 between 0 and 1
     2699      CoinBigIndex start[5] = { 0, 4, 7, 8, 9 };
     2700      int length[5] = { 2, 3, 1, 1, 1 };
     2701      int rows[11] = { 0, 2, -1, -1, 0, 1, 2, 0, 1, 2 };
     2702      double elements[11] = { 7.0, 2.0, 1.0e10, 1.0e10, -2.0, 1.0, -2.0, 1, 1, 1 };
     2703      CoinPackedMatrix matrix(true, 3, 5, 8, elements, rows, start, length);
     2704      // by row
     2705      matrix.reverseOrdering();
     2706      const double *element = matrix.getElements();
     2707      const int *column = matrix.getIndices();
     2708      const CoinBigIndex *rowStart = matrix.getVectorStarts();
     2709      const int *rowLength = matrix.getVectorLengths();
     2710
     2711      // rim data
     2712      //double objective[5]={-4.0,1.0,0.0,0.0,0.0};
     2713      double rowLower[3] = { 14.0, 3.0, 3.0 };
     2714      double rowUpper[3] = { 14.0, 3.0, 3.0 };
     2715      //double columnLower[5]={0.0,0.0,0.0,0.0,0.0};
     2716      //double columnUpper[5]={100.0,100.0,100.0,100.0,100.0};
     2717
     2718      for (int i = 0; i < 3; i++) {
     2719        sub.addRow(rowLength[i], column + rowStart[i],
     2720          element + rowStart[i], rowLower[i], rowUpper[i]);
     2721      }
     2722      //for (int i=0;i<5;i++) {
     2723      //sub.setColumnBounds(i,columnLower[i],columnUpper[i]);
     2724      //sub.setColumnObjective(i,objective[i]);
     2725      //}
     2726      sub.convertMatrix();
     2727    }
     2728    // Top
     2729    CoinModel top;
     2730    {
     2731      // matrix data
     2732      CoinBigIndex start[5] = { 0, 2, 4, 6, 8 };
     2733      int length[5] = { 2, 2, 2, 2, 2 };
     2734      int rows[10] = { 0, 1, 0, 1, 0, 1, 0, 1, 0, 1 };
     2735      double elements[10] = { 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 };
     2736      CoinPackedMatrix matrix(true, 2, 5, 8, elements, rows, start, length);
     2737      // by row
     2738      matrix.reverseOrdering();
     2739      const double *element = matrix.getElements();
     2740      const int *column = matrix.getIndices();
     2741      const CoinBigIndex *rowStart = matrix.getVectorStarts();
     2742      const int *rowLength = matrix.getVectorLengths();
     2743
     2744      // rim data
     2745      double objective[5] = { -4.0, 1.0, 0.0, 0.0, 0.0 };
     2746      //double rowLower[3]={14.0,3.0,3.0};
     2747      //double rowUpper[3]={14.0,3.0,3.0};
     2748      double columnLower[5] = { 0.0, 0.0, 0.0, 0.0, 0.0 };
     2749      double columnUpper[5] = { 100.0, 100.0, 100.0, 100.0, 100.0 };
     2750
     2751      for (int i = 0; i < 2; i++) {
     2752        top.addRow(rowLength[i], column + rowStart[i],
     2753          element + rowStart[i],
     2754          -COIN_DBL_MAX, COIN_DBL_MAX);
     2755      }
     2756      for (int i = 0; i < 5; i++) {
     2757        top.setColumnBounds(i, columnLower[i], columnUpper[i]);
     2758        top.setColumnObjective(i, objective[i]);
     2759      }
     2760      top.convertMatrix();
     2761    }
     2762    // Create a structured model
     2763    CoinStructuredModel structured;
     2764    int numberBlocks = 5;
     2765    for (int i = 0; i < numberBlocks; i++) {
     2766      std::string topName = "row_master";
     2767      std::string blockName = "block_";
     2768      char bName = static_cast< char >('a' + static_cast< char >(i));
     2769      blockName.append(1, bName);
     2770      structured.addBlock(topName, blockName, top);
     2771      structured.addBlock(blockName, blockName, sub);
     2772    }
     2773    // Set rhs on first block
     2774    CoinModel *first = structured.coinBlock(0);
     2775    for (int i = 0; i < 2; i++) {
     2776      first->setRowLower(i, 0.0);
     2777      first->setRowUpper(i, 100.0);
     2778    }
     2779    // Refresh whats set
     2780    structured.refresh(0);
     2781    // Could perturb stuff, but for first go don't bother
     2782    ClpSimplex fullModel;
     2783    // There is no original stuff set - think
     2784    fullModel.loadProblem(structured, false);
     2785    fullModel.dual();
     2786    fullModel.dropNames();
     2787    fullModel.writeMps("test.mps");
     2788    // Make up very simple nested model - not realistic
     2789    // Create a structured model
     2790    CoinStructuredModel structured2;
     2791    numberBlocks = 3;
     2792    for (int i = 0; i < numberBlocks; i++) {
     2793      std::string blockName = "block_";
     2794      char bName = static_cast< char >('a' + static_cast< char >(i));
     2795      blockName.append(1, bName);
     2796      structured2.addBlock(blockName, blockName, structured);
     2797    }
     2798    fullModel.loadProblem(structured2, false);
     2799    fullModel.dual();
     2800    fullModel.dropNames();
     2801    fullModel.writeMps("test2.mps");
     2802  }
    28062803#endif
    28072804}
     2805
     2806/* vi: softtabstop=2 shiftwidth=2 expandtab tabstop=2
     2807*/
Note: See TracChangeset for help on using the changeset viewer.