Ignore:
Timestamp:
Jun 19, 2011 1:23:14 PM (8 years ago)
Author:
stefan
Message:

sync with trunk rev1674

Location:
stable/2.7
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • stable/2.7

    • Property svn:externals
      •  

        old new  
        1 BuildTools  https://projects.coin-or.org/svn/BuildTools/stable/0.7
        2 ThirdParty/ASL  https://projects.coin-or.org/svn/BuildTools/ThirdParty/ASL/stable/1.2
        3 ThirdParty/Blas  https://projects.coin-or.org/svn/BuildTools/ThirdParty/Blas/stable/1.2
        4 ThirdParty/Lapack  https://projects.coin-or.org/svn/BuildTools/ThirdParty/Lapack/stable/1.3
        5 ThirdParty/Glpk  https://projects.coin-or.org/svn/BuildTools/ThirdParty/Glpk/stable/1.8
        6 Data/Sample  https://projects.coin-or.org/svn/Data/Sample/stable/1.2
        7 CoinUtils  https://projects.coin-or.org/svn/CoinUtils/stable/2.8/CoinUtils
        8 Cgl  https://projects.coin-or.org/svn/Cgl/stable/0.57/Cgl
        9 Clp  https://projects.coin-or.org/svn/Clp/stable/1.14/Clp
        10 Osi  https://projects.coin-or.org/svn/Osi/stable/0.105/Osi
         1BuildTools        https://projects.coin-or.org/svn/BuildTools/stable/0.7
         2ThirdParty/ASL    https://projects.coin-or.org/svn/BuildTools/ThirdParty/ASL/stable/1.2
         3ThirdParty/Blas   https://projects.coin-or.org/svn/BuildTools/ThirdParty/Blas/stable/1.2
         4ThirdParty/Lapack https://projects.coin-or.org/svn/BuildTools/ThirdParty/Lapack/stable/1.3
         5ThirdParty/Glpk   https://projects.coin-or.org/svn/BuildTools/ThirdParty/Glpk/stable/1.8
         6ThirdParty/Metis  https://projects.coin-or.org/svn/BuildTools/ThirdParty/Metis/stable/1.2
         7ThirdParty/Mumps  https://projects.coin-or.org/svn/BuildTools/ThirdParty/Mumps/stable/1.4
         8Data/Sample       https://projects.coin-or.org/svn/Data/Sample/stable/1.2
         9Data/miplib3      https://projects.coin-or.org/svn/Data/miplib3/stable/1.2
         10CoinUtils         https://projects.coin-or.org/svn/CoinUtils/stable/2.8/CoinUtils
         11Cgl               https://projects.coin-or.org/svn/Cgl/stable/0.57/Cgl
         12Clp               https://projects.coin-or.org/svn/Clp/stable/1.14/Clp
         13Osi               https://projects.coin-or.org/svn/Osi/stable/0.105/Osi
    • Property svn:mergeinfo changed
      /trunk (added)merged: 1578,​1582-1587,​1589-1600,​1603-1614,​1620-1626,​1631-1633,​1635-1636,​1638-1646,​1650-1652,​1654-1658,​1660-1663,​1665-1671,​1673-1674
  • stable/2.7/Cbc

  • stable/2.7/Cbc/src/unitTestClp.cpp

    r1573 r1675  
    88#include <string>
    99#include <iostream>
     10#include <iomanip>
    1011
    1112#include "CoinTime.hpp"
     13#include "CoinFileIO.hpp"
    1214#include "CbcModel.hpp"
    1315#include "CbcHeuristic.hpp"
     
    6163        return true;
    6264    }
    63 #ifdef COIN_HAS_ZLIB
    64     if (CbcTestFile(fname + ".gz")) {
    65         return true;
    66     }
    67     if (CbcTestFile(fname + ".mps.gz")) {
     65    if (CoinFileInput::haveGzipSupport()) {
     66      if (CbcTestFile(fname + ".gz")) {
     67        return true;
     68      }
     69      if (CbcTestFile(fname + ".mps.gz")) {
    6870        fname += ".mps";
    6971        return true;
    70     }
    71     if (CbcTestFile(fname + ".MPS.gz")) {
     72      }
     73      if (CbcTestFile(fname + ".MPS.gz")) {
    7274        fname += ".MPS";
    7375        return true;
    74     }
    75     if (CbcTestFile(fname + ".MPS.GZ")) {
     76      }
     77      if (CbcTestFile(fname + ".MPS.GZ")) {
    7678        fname += ".MPS";
    7779        return true;
    78     }
    79 #endif
     80      }
     81    }
     82    if (CoinFileInput::haveBzip2Support()) {
     83      if (CbcTestFile(fname + ".bz2")) {
     84        return true;
     85      }
     86      if (CbcTestFile(fname + ".mps.bz2")) {
     87        fname += ".mps";
     88        return true;
     89      }
     90      if (CbcTestFile(fname + ".MPS.bz2")) {
     91        fname += ".MPS";
     92        return true;
     93      }
     94      if (CbcTestFile(fname + ".MPS.BZ2")) {
     95        fname += ".MPS";
     96        return true;
     97      }
     98    }
    8099    return false;
    81100}
    82101
    83102//#############################################################################
    84 // testSwitch -2 unitTest, -1 normal (==2)
    85 int CbcClpUnitTest (const CbcModel & saveModel, std::string& dirMiplib,
    86                     int testSwitch,
    87                     double * stuff)
     103/*
     104  jjf: testSwitch -2 unitTest, -1 normal (==2)
     105
     106  MiplibTest might be more appropriate.
     107
     108  TestSwitch and stuff[6] together control how much of miplib is executed:
     109    For testSwitch set to:
     110      -2: solve p0033 and p0201 only (the unit test)
     111      -1: solve miplib sets #0 and #1
     112       0: solve nothing
     113       k: execute sets j:k, where j is determined by the value of stuff[6]
     114  The last parameter of PUSH_MPS specifies the test set membership.
     115
     116  For -miplib, -extra2 sets testSwitch, -extra3 sets stuff[6]. The command
     117    cbc -extra2 -2 -miplib
     118  will execute the unit test on the miplib directory.
     119
     120  dirMiplib should end in the directory separator character for the platform.
     121
     122  If you want to activate the row cut debugger for a given problem, change the
     123  last parameter of the PUSH_MPS macro for the problem to true.
     124
     125  Returns 0 if all goes well, -1 if the Miplib directory is missing, otherwise
     126     100*(number with bad objective)+(number that exceeded node limit)
     127*/
     128int CbcClpUnitTest (const CbcModel &saveModel, const std::string &dirMiplib,
     129                    int testSwitch, const double *stuff)
    88130{
    89     // Stop Windows popup
    90     WindowsErrorPopupBlocker();
    91     unsigned int m ;
    92 
    93     // Set directory containing miplib data files.
    94     std::string test1 = dirMiplib + "p0033";
    95     // See if files exist
    96     bool doTest = CbcTestMpsFile(test1);
    97 
    98     if (!doTest) {
    99         printf("Not doing miplib run as can't find mps files - ? .gz without libz\n");
    100         return -1;
    101     }
    102     /*
    103       Vectors to hold test problem names and characteristics. The objective value
    104       after optimization (objValue) must agree to the specified tolerance
    105       (objValueTol).
    106     */
     131  // Stop Windows popup
     132  WindowsErrorPopupBlocker() ;
     133  unsigned int m ;
     134
     135  // Do an existence check.
     136  std::string test1 = dirMiplib+"p0033";
     137  bool doTest = CbcTestMpsFile(test1);
     138  if (!doTest) {
     139    std::cout
     140      << "Not doing miplib run as can't find mps files." << std::endl
     141      << "Perhaps you're trying to read gzipped (.gz) files without libz?"
     142      << std::endl ;
     143    return (-1) ;
     144  }
     145  int dfltPrecision = static_cast<int>(std::cout.precision()) ;
     146/*
     147  Set the range of problems to be tested. testSwitch = -2 is special and is
     148  picked up below.
     149*/
     150    int loSet = 0 ;
     151    int hiSet = 0 ;
     152    if (testSwitch == -1) {
     153      loSet = 0 ;
     154      hiSet = 1 ;
     155    } else if (testSwitch >= 0) {
     156      loSet = static_cast<int>(stuff[6]) ;
     157      hiSet = testSwitch ;
     158      std::cout
     159        << "Solving miplib problems in sets " << loSet
     160        << ":" << hiSet << "." << std::endl ;
     161    }
     162/*
     163  Vectors to hold test problem names and characteristics.
     164*/
    107165    std::vector<std::string> mpsName ;
    108166    std::vector<int> nRows ;
     
    111169    std::vector<double> objValue ;
    112170    std::vector<int> testSet ;
    113     /*
    114       And a macro to make the vector creation marginally readable.
    115     */
     171    std::vector<bool> rowCutDebugger ;
     172/*
     173  A macro to make the vector creation marginally readable. Parameters are
     174  name, rows, columns, integer objective, continuous objective, set ID,
     175  row cut debugger
     176
     177  To enable the row cut debugger for a given problem, change the last
     178  parameter to true. Don't forget to turn it off before committing changes!
     179*/
    116180#define PUSH_MPS(zz_mpsName_zz,\
    117181                 zz_nRows_zz,zz_nCols_zz,zz_objValue_zz,zz_objValueC_zz, \
    118                  zz_testSet_zz) \
     182                 zz_testSet_zz, zz_rcDbg_zz) \
    119183  mpsName.push_back(zz_mpsName_zz) ; \
    120184  nRows.push_back(zz_nRows_zz) ; \
     
    122186  objValueC.push_back(zz_objValueC_zz) ; \
    123187  testSet.push_back(zz_testSet_zz) ; \
    124   objValue.push_back(zz_objValue_zz) ;
    125     int loSwitch = 0;
    126     if (testSwitch == -2) {
    127         PUSH_MPS("p0033", 16, 33, 3089, 2520.57, 0);
    128         PUSH_MPS("p0201", 133, 201, 7615, 6875.0, 0);
    129         testSwitch = 0;
     188  objValue.push_back(zz_objValue_zz) ; \
     189  rowCutDebugger.push_back(zz_rcDbg_zz) ;
     190/*
     191  Push the miplib problems. Except for -2 (unitTest), push all, even if we're
     192  not going to do all of them.
     193*/
     194  if (testSwitch == -2) {
     195      PUSH_MPS("p0033", 16, 33, 3089, 2520.57, 0, false);
     196      PUSH_MPS("p0201", 133, 201, 7615, 6875.0, 0, false);
     197      // PUSH_MPS("flugpl", 18, 18, 1201500, 1167185.7, 0, false);
     198  } else {
     199/*
     200  Load up the problem vector. Note that the row counts here include the
     201  objective function.
     202*/
     203    PUSH_MPS("10teams", 230, 2025, 924, 917, 1, false);
     204    PUSH_MPS("air03", 124, 10757, 340160, 338864.25, 0, false);
     205    PUSH_MPS("air04", 823, 8904, 56137, 55535.436, 2, false);
     206    PUSH_MPS("air05", 426, 7195, 26374, 25877.609, 2, false);
     207    PUSH_MPS("arki001", 1048, 1388, 7580813.0459, 7579599.80787, 7, false);
     208    PUSH_MPS("bell3a", 123, 133, 878430.32, 862578.64, 0, false);
     209    PUSH_MPS("bell5", 91, 104, 8966406.49, 8608417.95, 1, false);
     210    PUSH_MPS("blend2", 274, 353, 7.598985, 6.9156751140, 0, false);
     211    PUSH_MPS("cap6000", 2176, 6000, -2451377, -2451537.325, 1, false);
     212    PUSH_MPS("dano3mip", 3202, 13873, 728.1111, 576.23162474, 7, false);
     213    PUSH_MPS("danoint", 664, 521, 65.67, 62.637280418, 6, false);
     214    PUSH_MPS("dcmulti", 290, 548, 188182, 183975.5397, 0, false);
     215    PUSH_MPS("dsbmip", 1182, 1886, -305.19817501, -305.19817501, 0, false);
     216    PUSH_MPS("egout", 98, 141, 568.101, 149.589, 0, false);
     217    PUSH_MPS("enigma", 21, 100, 0.0, 0.0, 0, false);
     218    PUSH_MPS("fast0507", 507, 63009, 174, 172.14556668, 5, false);
     219    PUSH_MPS("fiber", 363, 1298, 405935.18000, 156082.51759, 0, false);
     220    PUSH_MPS("fixnet6", 478, 878, 3983, 1200.88, 1, false);
     221    PUSH_MPS("flugpl", 18, 18, 1201500, 1167185.7, 0, false);
     222    PUSH_MPS("gen", 780, 870, 112313, 112130.0, 0, false);
     223    PUSH_MPS("gesa2", 1392, 1224, 25779856.372, 25476489.678, 1, false);
     224    PUSH_MPS("gesa2_o", 1248, 1224, 25779856.372, 25476489.678, 1, false);
     225    PUSH_MPS("gesa3", 1368, 1152, 27991042.648, 27833632.451, 0, false);
     226    PUSH_MPS("gesa3_o", 1224, 1152, 27991042.648, 27833632.451, 0, false);
     227    PUSH_MPS("gt2", 29, 188, 21166.000, 13460.233074, 0, false);
     228    PUSH_MPS("harp2", 112, 2993, -73899798.00, -74353341.502, 6, false);
     229    PUSH_MPS("khb05250", 101, 1350, 106940226, 95919464.0, 0, false);
     230    PUSH_MPS("l152lav", 97, 1989, 4722, 4656.36, 1, false);
     231    PUSH_MPS("lseu", 28, 89, 1120, 834.68, 0, false);
     232    PUSH_MPS("mas74", 13, 151, 11801.18573, 10482.79528, 3, false);
     233    PUSH_MPS("mas76", 12, 151, 40005.05414, 38893.9036, 2, false);
     234    PUSH_MPS("misc03", 96, 160, 3360, 1910., 0, false);
     235    PUSH_MPS("misc06", 820, 1808, 12850.8607, 12841.6, 0, false);
     236    PUSH_MPS("misc07", 212, 260, 2810, 1415.0, 1, false);
     237    PUSH_MPS("mitre", 2054, 10724, 115155, 114740.5184, 1, false);
     238    PUSH_MPS("mkc", 3411, 5325, -553.75, -611.85, 7, false); // suboptimal
     239    PUSH_MPS("mod008", 6, 319, 307, 290.9, 0, false);
     240    PUSH_MPS("mod010", 146, 2655, 6548, 6532.08, 0, false);
     241    PUSH_MPS("mod011", 4480, 10958, -54558535, -62121982.55, 2, false);
     242    PUSH_MPS("modglob", 291, 422, 20740508, 20430947., 2, false);
     243    PUSH_MPS("noswot", 182, 128, -43, -43.0, 6, false);
     244    PUSH_MPS("nw04", 36, 87482, 16862, 16310.66667, 1, false);
     245    PUSH_MPS("p0033", 16, 33, 3089, 2520.57, 0, false);
     246    PUSH_MPS("p0201", 133, 201, 7615, 6875.0, 0, false);
     247    PUSH_MPS("p0282", 241, 282, 258411, 176867.50, 0, false);
     248    PUSH_MPS("p0548", 176, 548, 8691, 315.29, 0, false);
     249    PUSH_MPS("p2756", 755, 2756, 3124, 2688.75, 0, false);
     250    PUSH_MPS("pk1", 45, 86, 11.0, 0.0, 2, false);
     251    PUSH_MPS("pp08a", 136, 240, 7350.0, 2748.3452381, 1, false);
     252    PUSH_MPS("pp08aCUTS", 246, 240, 7350.0, 5480.6061563, 1, false);
     253    PUSH_MPS("qiu", 1192, 840, -132.873137, -931.638857, 3, false);
     254    PUSH_MPS("qnet1", 503, 1541, 16029.692681, 14274.102667, 0, false);
     255    PUSH_MPS("qnet1_o", 456, 1541, 16029.692681, 12095.571667, 0, false);
     256    PUSH_MPS("rentacar", 6803, 9557, 30356761, 28806137.644, 0, false);
     257    PUSH_MPS("rgn", 24, 180, 82.1999, 48.7999, 0, false);
     258    PUSH_MPS("rout", 291, 556, 1077.56, 981.86428571, 3, false);
     259    PUSH_MPS("set1ch", 492, 712, 54537.75, 32007.73, 5, false);
     260    PUSH_MPS("seymour", 4944, 1372, 423, 403.84647413, 7, false);
     261    PUSH_MPS("seymour_1", 4944, 1372, 410.76370, 403.84647413, 5, false);
     262    PUSH_MPS("stein27", 118, 27, 18, 13.0, 0, false);
     263    PUSH_MPS("stein45", 331, 45, 30, 22.0, 1, false);
     264    PUSH_MPS("swath", 884, 6805, 497.603, 334.4968581, 7, false);
     265    PUSH_MPS("vpm1", 234, 378, 20, 15.4167, 0, false);
     266    PUSH_MPS("vpm2", 234, 378, 13.75, 9.8892645972, 0, false);
     267  }
     268#undef PUSH_MPS
     269
     270/*
     271  Normally the problems are executed in order. Define RANDOM_ORDER below to
     272  randomize.
     273
     274  #define RANDOM_ORDER
     275*/
     276  int which[100];
     277  int nLoop = static_cast<int>(mpsName.size());
     278  assert (nLoop <= 100);
     279  for (int i = 0; i < nLoop; i++) which[i] = i;
     280
     281# ifdef RANDOM_ORDER
     282  unsigned int iTime = static_cast<unsigned int>(CoinGetTimeOfDay()-1.256e9);
     283  std::cout << "Time (seed) " << iTime << "." << std::endl ;
     284  double sort[100];
     285  CoinDrand48(true,iTime);
     286  for (int i = 0; i < nLoop; i++) sort[i] = CoinDrand48();
     287  CoinSort_2(sort,sort+nLoop,which);
     288# endif
     289
     290  int problemCnt = 0;
     291  for (m = 0 ; m < mpsName.size() ; m++) {
     292    int setID = testSet[m];
     293    if (loSet <= setID && setID <= hiSet) problemCnt++;
     294  }
     295
     296  int numberFailures = 0;
     297  int numberAttempts = 0;
     298  int numProbSolved = 0;
     299  double timeTaken = 0.0;
     300
     301//#define CLP_FACTORIZATION_INSTRUMENT
     302# ifdef CLP_FACTORIZATION_INSTRUMENT
     303  double timeTakenFac = 0.0;
     304# endif
     305
     306/*
     307  Open the main loop to step through the MPS problems.
     308*/
     309  for (unsigned int mw = 0 ; mw < mpsName.size() ; mw++) {
     310    m = which[mw];
     311    int setID = testSet[m];
     312    // Skip if problem is not in specified problem set(s)
     313    if (!(loSet <= setID && setID <= hiSet)) continue ;
     314   
     315    numberAttempts++;
     316    std::cout << "  processing mps file: " << mpsName[m]
     317              << " (" << numberAttempts << " out of "
     318              << problemCnt << ")" << std::endl ;
     319/*
     320  Stage 1: Read the MPS and make sure the size of the constraint matrix
     321           is correct.
     322*/
     323    CbcModel *model = new CbcModel(saveModel) ;
     324
     325    std::string fn = dirMiplib+mpsName[m] ;
     326    if (!CbcTestMpsFile(fn)) {
     327      std::cout << "ERROR: Cannot find MPS file " << fn << "." << std::endl ;
     328      continue;
     329    }
     330    model->solver()->readMps(fn.c_str(), "") ;
     331    assert(model->getNumRows() == nRows[m]) ;
     332    assert(model->getNumCols() == nCols[m]) ;
     333
     334    // Careful! We're initialising for the benefit of other code.
     335    CoinDrand48(true,1234567);
     336    //printf("RAND1 %g %g\n",CoinDrand48(true,1234567),model->randomNumberGenerator()->randomDouble());
     337    //printf("RAND1 %g\n",CoinDrand48(true,1234567));
     338
     339    // Higher limits for the serious problems.
     340    int testMaximumNodes = 200000;
     341    if (hiSet > 1)
     342        testMaximumNodes = 20000000;
     343    if (model->getMaximumNodes() > testMaximumNodes) {
     344        model->setMaximumNodes(testMaximumNodes);
     345    }
     346/*
     347  Stage 2: Call solver to solve the problem.
     348*/
     349
     350#   ifdef CLP_FACTORIZATION_INSTRUMENT
     351    extern double factorization_instrument(int type);
     352    double facTime1 = factorization_instrument(0);
     353    std::cout
     354      << "Factorization - initial solve " << facTime1 << " seconds."
     355      << std::endl ;
     356      timeTakenFac += facTime1;
     357#   endif
     358
     359    double startTime = CoinCpuTime()+CoinCpuTimeJustChildren();
     360
     361    // Setup specific to clp
     362    OsiClpSolverInterface *siClp =
     363        dynamic_cast<OsiClpSolverInterface *>(model->solver()) ;
     364    ClpSimplex *modelC = NULL;
     365    if (siClp) {
     366      modelC = siClp->getModelPtr();
     367      ClpMatrixBase * matrix = modelC->clpMatrix();
     368      ClpPackedMatrix * clpMatrix = dynamic_cast< ClpPackedMatrix*>(matrix);
     369      if (stuff && stuff[9] && clpMatrix) {
     370        // vector matrix!
     371        clpMatrix->makeSpecialColumnCopy();
     372      }
     373
     374#     ifdef JJF_ZERO
     375      if (clpMatrix) {
     376        int numberRows = clpMatrix->getNumRows();
     377        int numberColumns = clpMatrix->getNumCols();
     378        double * elements = clpMatrix->getMutableElements();
     379        const int * row = clpMatrix->getIndices();
     380        const CoinBigIndex * columnStart = clpMatrix->getVectorStarts();
     381        const int * columnLength = clpMatrix->getVectorLengths();
     382        double * smallest = new double [numberRows];
     383        double * largest = new double [numberRows];
     384        char * flag = new char [numberRows];
     385        CoinZeroN(flag, numberRows);
     386        for (int i = 0; i < numberRows; i++) {
     387            smallest[i] = COIN_DBL_MAX;
     388            largest[i] = 0.0;
     389        }
     390        for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
     391            bool isInteger = modelC->isInteger(iColumn);
     392            CoinBigIndex j;
     393            for (j = columnStart[iColumn];
     394                    j < columnStart[iColumn] + columnLength[iColumn]; j++) {
     395                int iRow = row[j];
     396                double value = fabs(elements[j]);
     397                if (!isInteger)
     398                    flag[iRow] = 1;
     399                smallest[iRow] = CoinMin(smallest[iRow], value);
     400                largest[iRow] = CoinMax(largest[iRow], value);
     401            }
     402        }
     403        double * rowLower = modelC->rowLower();
     404        double * rowUpper = modelC->rowUpper();
     405        bool changed = false;
     406        for (int i = 0; i < numberRows; i++) {
     407            if (flag[i] && smallest[i] > 10.0 && false) {
     408                smallest[i] = 1.0 / smallest[i];
     409                if (rowLower[i] > -1.0e20)
     410                    rowLower[i] *= smallest[i];
     411                if (rowUpper[i] < 1.0e20)
     412                    rowUpper[i] *= smallest[i];
     413                changed = true;
     414            } else {
     415                smallest[i] = 0.0;
     416            }
     417        }
     418        if (changed) {
     419            printf("SCALED\n");
     420            for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
     421                CoinBigIndex j;
     422                for (j = columnStart[iColumn];
     423                        j < columnStart[iColumn] + columnLength[iColumn]; j++) {
     424                    int iRow = row[j];
     425                    if (smallest[iRow])
     426                        elements[j] *= smallest[iRow];
     427                }
     428            }
     429        }
     430        delete [] smallest;
     431        delete [] largest;
     432        delete [] flag;
     433      }
     434#     endif    // JJF_ZERO
     435
     436      model->checkModel();
     437      modelC->tightenPrimalBounds(0.0, 0, true);
     438      model->initialSolve();
     439      if (modelC->dualBound() == 1.0e10) {
     440        // user did not set - so modify
     441        // get largest scaled away from bound
     442        ClpSimplex temp = *modelC;
     443        temp.dual(0, 7);
     444        double largestScaled = 1.0e-12;
     445        double largest = 1.0e-12;
     446        int numberRows = temp.numberRows();
     447        const double * rowPrimal = temp.primalRowSolution();
     448        const double * rowLower = temp.rowLower();
     449        const double * rowUpper = temp.rowUpper();
     450        const double * rowScale = temp.rowScale();
     451        int iRow;
     452        for (iRow = 0; iRow < numberRows; iRow++) {
     453          double value = rowPrimal[iRow];
     454          double above = value - rowLower[iRow];
     455          double below = rowUpper[iRow] - value;
     456          if (above < 1.0e12) {
     457              largest = CoinMax(largest, above);
     458          }
     459          if (below < 1.0e12) {
     460              largest = CoinMax(largest, below);
     461          }
     462          if (rowScale) {
     463              double multiplier = rowScale[iRow];
     464              above *= multiplier;
     465              below *= multiplier;
     466          }
     467          if (above < 1.0e12) {
     468              largestScaled = CoinMax(largestScaled, above);
     469          }
     470          if (below < 1.0e12) {
     471              largestScaled = CoinMax(largestScaled, below);
     472          }
     473        }
     474
     475        int numberColumns = temp.numberColumns();
     476        const double * columnPrimal = temp.primalColumnSolution();
     477        const double * columnLower = temp.columnLower();
     478        const double * columnUpper = temp.columnUpper();
     479        const double * columnScale = temp.columnScale();
     480        int iColumn;
     481        for (iColumn = 0; iColumn < numberColumns; iColumn++) {
     482          double value = columnPrimal[iColumn];
     483          double above = value - columnLower[iColumn];
     484          double below = columnUpper[iColumn] - value;
     485          if (above < 1.0e12) {
     486              largest = CoinMax(largest, above);
     487          }
     488          if (below < 1.0e12) {
     489              largest = CoinMax(largest, below);
     490          }
     491          if (columnScale) {
     492              double multiplier = 1.0 / columnScale[iColumn];
     493              above *= multiplier;
     494              below *= multiplier;
     495          }
     496          if (above < 1.0e12) {
     497              largestScaled = CoinMax(largestScaled, above);
     498          }
     499          if (below < 1.0e12) {
     500              largestScaled = CoinMax(largestScaled, below);
     501          }
     502        }
     503        std::cout << "Largest (scaled) away from bound " << largestScaled
     504                  << " unscaled " << largest << std::endl;
     505#       ifdef JJF_ZERO
     506        modelC->setDualBound(CoinMax(1.0001e8,
     507                                 CoinMin(1000.0*largestScaled,1.00001e10)));
     508#       else
     509        modelC->setDualBound(CoinMax(1.0001e9,
     510                                 CoinMin(1000.0*largestScaled,1.0001e10)));
     511#       endif
     512      }
     513    }    // end clp-specific setup
     514/*
     515  Cut passes: For small models (n < 500) always do 100 passes, if possible
     516  (-100). For larger models, use minimum drop to stop (100, 20).
     517*/
     518    model->setMinimumDrop(CoinMin(5.0e-2,
     519                      fabs(model->getMinimizationObjValue())*1.0e-3+1.0e-4));
     520    if (CoinAbs(model->getMaximumCutPassesAtRoot()) <= 100) {
     521        if (model->getNumCols() < 500) {
     522          model->setMaximumCutPassesAtRoot(-100);
     523        } else if (model->getNumCols() < 5000) {
     524          model->setMaximumCutPassesAtRoot(100);
     525        } else {
     526          model->setMaximumCutPassesAtRoot(20);
     527        }
     528    }
     529    // If defaults then increase trust for small models
     530    if (model->numberStrong() == 5 && model->numberBeforeTrust() == 10) {
     531        int numberColumns = model->getNumCols();
     532        if (numberColumns <= 50) {
     533            model->setNumberBeforeTrust(1000);
     534        } else if (numberColumns <= 100) {
     535            model->setNumberBeforeTrust(100);
     536        } else if (numberColumns <= 300) {
     537            model->setNumberBeforeTrust(50);
     538        }
     539    }
     540    //if (model->getNumCols()>=500) {
     541    // switch off Clp stuff
     542    //model->setFastNodeDepth(-1);
     543    //}
     544/*
     545  Activate the row cut debugger, if requested.
     546*/
     547    if (rowCutDebugger[m] == true) {
     548      std::string probName ;
     549      model->solver()->getStrParam(OsiProbName,probName) ;
     550      model->solver()->activateRowCutDebugger(probName.c_str()) ;
     551      if (model->solver()->getRowCutDebugger())
     552        std::cout << "Row cut debugger activated for " ;
     553      else
     554        std::cout << "Failed to activate row cut debugger for " ;
     555      std::cout << mpsName[m] << "." << std::endl ;
     556    }
     557    setCutAndHeuristicOptions(*model) ;
     558/*
     559  More clp-specific setup.
     560*/
     561    if (siClp) {
     562#     ifdef CLP_MULTIPLE_FACTORIZATIONS
     563      if (!modelC->factorization()->isDenseOrSmall()) {
     564          int denseCode = stuff ? static_cast<int> (stuff[4]) : -1;
     565          int smallCode = stuff ? static_cast<int> (stuff[10]) : -1;
     566          if (stuff && stuff[8] >= 1) {
     567              if (denseCode < 0)
     568                  denseCode = 40;
     569              if (smallCode < 0)
     570                  smallCode = 40;
     571          }
     572          if (denseCode > 0)
     573              modelC->factorization()->setGoDenseThreshold(denseCode);
     574          if (smallCode > 0)
     575              modelC->factorization()->setGoSmallThreshold(smallCode);
     576          if (denseCode >= modelC->numberRows()) {
     577              //printf("problem going dense\n");
     578              //modelC->factorization()->goDenseOrSmall(modelC->numberRows());
     579          }
     580      }
     581#     endif
     582      if (stuff && stuff[8] >= 1) {
     583          if (modelC->numberColumns() + modelC->numberRows() <= 10000 &&
     584                  model->fastNodeDepth() == -1)
     585              model->setFastNodeDepth(-9);
     586      }
     587    }
     588    //OsiObject * obj = new CbcBranchToFixLots(model,0.3,0.0,3,3000003);
     589    //model->addObjects(1,&obj);
     590    //delete obj;
     591/*
     592  Finally, the actual call to solve the MIP with branch-and-cut.
     593*/
     594    model->branchAndBound();
     595
     596#   ifdef CLP_FACTORIZATION_INSTRUMENT
     597    double facTime = factorization_instrument(0);
     598    std::cout << "Factorization " << facTime << " seconds." << std::endl ,
     599    timeTakenFac += facTime;
     600#   endif
     601
     602/*
     603  Stage 3: Do the statistics and check the answer.
     604*/
     605    double timeOfSolution = CoinCpuTime()+CoinCpuTimeJustChildren()-startTime;
     606    std::cout
     607      << "Cuts at root node changed objective from "
     608      << model->getContinuousObjective() << " to "
     609      << model->rootObjectiveAfterCuts() << std::endl ;
     610    int numberGenerators = model->numberCutGenerators();
     611    for (int iGenerator = 0 ; iGenerator < numberGenerators ; iGenerator++) {
     612      CbcCutGenerator *generator = model->cutGenerator(iGenerator);
     613#     ifdef CLIQUE_ANALYSIS
     614#     ifndef CLP_INVESTIGATE
     615      CglImplication *implication =
     616            dynamic_cast<CglImplication*>(generator->generator());
     617      if (implication) continue;
     618#     endif
     619#     endif
     620      std::cout
     621        << generator->cutGeneratorName() << " was tried "
     622        << generator->numberTimesEntered() << " times and created "
     623        << generator->numberCutsInTotal() << " cuts of which "
     624        << generator->numberCutsActive()
     625        << " were active after adding rounds of cuts";
     626      if (generator->timing())
     627        std::cout << " (" << generator->timeInCutGenerator() << " seconds)" ;
     628      std::cout << "." << std::endl;
     629    }
     630    std::cout
     631      << model->getNumberHeuristicSolutions()
     632      << " solutions found by heuristics." << std::endl ;
     633    int numberHeuristics = model->numberHeuristics();
     634    for (int iHeuristic = 0 ; iHeuristic < numberHeuristics ; iHeuristic++) {
     635      CbcHeuristic *heuristic = model->heuristic(iHeuristic);
     636      if (heuristic->numRuns()) {
     637        std::cout
     638          << heuristic->heuristicName() << " was tried "
     639          << heuristic->numRuns() << " times out of "
     640          << heuristic->numCouldRun() << " and created "
     641          << heuristic->numberSolutionsFound() << " solutions." << std::endl ;
     642        }
     643    }
     644/*
     645  Check for the correct answer.
     646*/
     647    if (!model->status()) {
     648
     649      double objActual = model->getObjValue() ;
     650      double objExpect = objValue[m] ;
     651      double tolerance = CoinMin(fabs(objActual),fabs(objExpect)) ;
     652      tolerance = CoinMax(1.0e-5,1.0e-5*tolerance) ;
     653      //CoinRelFltEq eq(1.0e-3) ;
     654
     655      std::cout
     656        << "cbc_clp (" << mpsName[m] << ") "
     657        << std::setprecision(10) << objActual ;
     658      if (fabs(objActual-objExpect) < tolerance) {
     659        std::cout << std::setprecision(dfltPrecision) << "; okay" ;
     660          numProbSolved++;
     661      } else  {
     662        std::cout
     663          << " != " << objExpect << std::setprecision(dfltPrecision)
     664          << "; error = " << fabs(objExpect-objActual) ;
     665          numberFailures++;
     666        //#ifdef COIN_DEVELOP
     667        //abort();
     668        //#endif
     669      }
    130670    } else {
    131         if (testSwitch == -1) {
    132             testSwitch = 1;
    133         } else {
    134             loSwitch = static_cast<int>(stuff[6]);
    135             printf("Solving miplib problems in sets >= %d and <=%d\n",
    136                    loSwitch, testSwitch);
    137         }
    138         /*
    139           Load up the problem vector. Note that the row counts here include the
    140           objective function.
    141         */
    142         // 0 for no test, 1 for some, 2 for many, 3 for all
    143         //PUSH_MPS("blend2",274,353,7.598985,6.9156751140,0);
    144         //PUSH_MPS("p2756",755,2756,3124,2688.75,0);
    145         //PUSH_MPS("seymour_1",4944,1372,410.7637014,404.35152,0);
    146         //PUSH_MPS("enigma",21,100,0.0,0.0,0);
    147         //PUSH_MPS("misc03",96,160,3360,1910.,0);
    148         //PUSH_MPS("p0201",133,201,7615,6875.0,0);
    149 #define HOWMANY 6
    150 #if HOWMANY
    151         PUSH_MPS("10teams", 230, 2025, 924, 917, 1);
    152         PUSH_MPS("air03", 124, 10757, 340160, 338864.25, 0);
    153         PUSH_MPS("air04", 823, 8904, 56137, 55535.436, 2);
    154         PUSH_MPS("air05", 426, 7195, 26374, 25877.609, 2);
    155         PUSH_MPS("arki001", 1048, 1388, 7580813.0459, 7579599.80787, 7);
    156         PUSH_MPS("bell3a", 123, 133, 878430.32, 862578.64, 0);
    157         PUSH_MPS("bell5", 91, 104, 8966406.49, 8608417.95, 1);
    158         PUSH_MPS("blend2", 274, 353, 7.598985, 6.9156751140, 0);
    159         PUSH_MPS("cap6000", 2176, 6000, -2451377, -2451537.325, 1);
    160         PUSH_MPS("dano3mip", 3202, 13873, 728.1111, 576.23162474, 7);
    161         PUSH_MPS("danoint", 664, 521, 65.67, 62.637280418, 6);
    162         PUSH_MPS("dcmulti", 290, 548, 188182, 183975.5397, 0);
    163         PUSH_MPS("dsbmip", 1182, 1886, -305.19817501, -305.19817501, 0);
    164         PUSH_MPS("egout", 98, 141, 568.101, 149.589, 0);
    165         PUSH_MPS("enigma", 21, 100, 0.0, 0.0, 0);
    166         PUSH_MPS("fast0507", 507, 63009, 174, 172.14556668, 5);
    167         PUSH_MPS("fiber", 363, 1298, 405935.18000, 156082.51759, 0);
    168         PUSH_MPS("fixnet6", 478, 878, 3983, 1200.88, 1);
    169         PUSH_MPS("flugpl", 18, 18, 1201500, 1167185.7, 0);
    170         PUSH_MPS("gen", 780, 870, 112313, 112130.0, 0);
    171         PUSH_MPS("gesa2", 1392, 1224, 25779856.372, 25476489.678, 1);
    172         PUSH_MPS("gesa2_o", 1248, 1224, 25779856.372, 25476489.678, 1);
    173         PUSH_MPS("gesa3", 1368, 1152, 27991042.648, 27833632.451, 0);
    174         PUSH_MPS("gesa3_o", 1224, 1152, 27991042.648, 27833632.451, 0);
    175         PUSH_MPS("gt2", 29, 188, 21166.000, 13460.233074, 0);
    176         PUSH_MPS("harp2", 112, 2993, -73899798.00, -74353341.502, 6);
    177         PUSH_MPS("khb05250", 101, 1350, 106940226, 95919464.0, 0);
    178         PUSH_MPS("l152lav", 97, 1989, 4722, 4656.36, 1);
    179         PUSH_MPS("lseu", 28, 89, 1120, 834.68, 0);
    180         PUSH_MPS("mas74", 13, 151, 11801.18573, 10482.79528, 3);
    181         PUSH_MPS("mas76", 12, 151, 40005.05414, 38893.9036, 2);
    182         PUSH_MPS("misc03", 96, 160, 3360, 1910., 0);
    183         PUSH_MPS("misc06", 820, 1808, 12850.8607, 12841.6, 0);
    184         PUSH_MPS("misc07", 212, 260, 2810, 1415.0, 1);
    185         PUSH_MPS("mitre", 2054, 10724, 115155, 114740.5184, 1);
    186         PUSH_MPS("mkc", 3411, 5325, -553.75, -611.85, 7); // this is suboptimal
    187         PUSH_MPS("mod008", 6, 319, 307, 290.9, 0);
    188         PUSH_MPS("mod010", 146, 2655, 6548, 6532.08, 0);
    189         PUSH_MPS("mod011", 4480, 10958, -54558535, -62121982.55, 2);
    190         PUSH_MPS("modglob", 291, 422, 20740508, 20430947., 2);
    191         PUSH_MPS("noswot", 182, 128, -43, -43.0, 6);
    192         PUSH_MPS("nw04", 36, 87482, 16862, 16310.66667, 1);
    193         PUSH_MPS("p0033", 16, 33, 3089, 2520.57, 0);
    194         PUSH_MPS("p0201", 133, 201, 7615, 6875.0, 0);
    195         PUSH_MPS("p0282", 241, 282, 258411, 176867.50, 0);
    196         PUSH_MPS("p0548", 176, 548, 8691, 315.29, 0);
    197         PUSH_MPS("p2756", 755, 2756, 3124, 2688.75, 0);
    198         PUSH_MPS("pk1", 45, 86, 11.0, 0.0, 2);
    199         PUSH_MPS("pp08a", 136, 240, 7350.0, 2748.3452381, 1);
    200         PUSH_MPS("pp08aCUTS", 246, 240, 7350.0, 5480.6061563, 1);
    201         PUSH_MPS("qiu", 1192, 840, -132.873137, -931.638857, 3);
    202         PUSH_MPS("qnet1", 503, 1541, 16029.692681, 14274.102667, 0);
    203         PUSH_MPS("qnet1_o", 456, 1541, 16029.692681, 12095.571667, 0);
    204         PUSH_MPS("rentacar", 6803, 9557, 30356761, 28806137.644, 0);
    205         PUSH_MPS("rgn", 24, 180, 82.1999, 48.7999, 0);
    206         PUSH_MPS("rout", 291, 556, 1077.56, 981.86428571, 3);
    207         PUSH_MPS("set1ch", 492, 712, 54537.75, 32007.73, 5);
    208         PUSH_MPS("seymour", 4944, 1372, 423, 403.84647413, 7);
    209         PUSH_MPS("seymour_1", 4944, 1372, 410.76370, 403.84647413, 5);
    210         PUSH_MPS("stein27", 118, 27, 18, 13.0, 0);
    211         PUSH_MPS("stein45", 331, 45, 30, 22.0, 1);
    212         PUSH_MPS("swath", 884, 6805, 497.603, 334.4968581, 7);
    213         PUSH_MPS("vpm1", 234, 378, 20, 15.4167, 0);
    214         PUSH_MPS("vpm2", 234, 378, 13.75, 9.8892645972, 0);
    215 #endif
    216     }
    217 #undef PUSH_MPS
    218 
    219     int numProbSolved = 0;
    220     double timeTaken = 0.0;
    221     //#define CLP_FACTORIZATION_INSTRUMENT
    222 #ifdef CLP_FACTORIZATION_INSTRUMENT
    223     double timeTakenFac = 0.0;
    224 #endif
    225     // Normally do in order
    226     int which[100];
    227     int nLoop = static_cast<int>(mpsName.size());
    228     assert (nLoop <= 100);
    229     for (int i = 0; i < nLoop; i++)
    230         which[i] = i;
    231     //#define RANDOM_ORDER
    232 #ifdef RANDOM_ORDER
    233     unsigned int iTime = static_cast<unsigned int>(CoinGetTimeOfDay() - 1.256e9);
    234     printf("Time %d\n", iTime);
    235     double sort[100];
    236     CoinDrand48(true, iTime);
    237     for (int i = 0; i < nLoop; i++)
    238         sort[i] = CoinDrand48();
    239     CoinSort_2(sort, sort + nLoop, which);
    240 #endif
    241     int numberFailures = 0;
    242     int numberAttempts = 0;
    243     int numberPossibleAttempts = 0;
    244     for (m = 0 ; m < mpsName.size() ; m++) {
    245         int test = testSet[m];
    246         if (testSwitch >= test && loSwitch <= test)
    247             numberPossibleAttempts++;
    248     }
    249 
    250     /*
    251       Open the main loop to step through the MPS problems.
    252     */
    253     for (unsigned int mw = 0 ; mw < mpsName.size() ; mw++) {
    254         m = which[mw];
    255         int test = testSet[m];
    256         if (testSwitch >= test && loSwitch <= test) {
    257             numberAttempts++;
    258             std::cout << "  processing mps file: " << mpsName[m]
    259                       << " (" << numberAttempts << " out of "
    260                       << numberPossibleAttempts << ")\n";
    261             /*
    262             Stage 1: Read the MPS
    263             and make sure the size of the constraint matrix is correct.
    264             */
    265             CbcModel * model = new CbcModel(saveModel);
    266 
    267             std::string fn = dirMiplib + mpsName[m] ;
    268             if (!CbcTestMpsFile(fn)) {
    269                 std::cout << "ERROR: Cannot find MPS file " << fn << "\n";
    270                 continue;
    271             }
    272             CoinDrand48(true, 1234567);
    273             //printf("RAND1 %g %g\n",CoinDrand48(true,1234567),model->randomNumberGenerator()->randomDouble());
    274             //printf("RAND1 %g\n",CoinDrand48(true,1234567));
    275             model->solver()->readMps(fn.c_str(), "") ;
    276             assert(model->getNumRows() == nRows[m]) ;
    277             assert(model->getNumCols() == nCols[m]) ;
    278 
    279             /*
    280               Stage 2: Call solver to solve the problem.  then check the return code and
    281                   objective.
    282             */
    283 
    284 #ifdef CLP_FACTORIZATION_INSTRUMENT
    285             extern double factorization_instrument(int type);
    286             double facTime1 = factorization_instrument(0);
    287             printf("Factorization - initial solve %g seconds\n",
    288                    facTime1);
    289             timeTakenFac += facTime1;
    290 #endif
    291             double startTime = CoinCpuTime() + CoinCpuTimeJustChildren();
    292             int testMaximumNodes = 200000;
    293             if (testSwitch > 1)
    294                 testMaximumNodes = 20000000;
    295             if (model->getMaximumNodes() > testMaximumNodes) {
    296                 model->setMaximumNodes(testMaximumNodes);
    297             }
    298             OsiClpSolverInterface * si =
    299                 dynamic_cast<OsiClpSolverInterface *>(model->solver()) ;
    300             ClpSimplex * modelC = NULL;
    301             if (si) {
    302                 // get clp itself
    303                 modelC = si->getModelPtr();
    304                 ClpMatrixBase * matrix = modelC->clpMatrix();
    305                 ClpPackedMatrix * clpMatrix = dynamic_cast< ClpPackedMatrix*>(matrix);
    306                 if (stuff && stuff[9] && clpMatrix) {
    307                     // vector matrix!
    308                     clpMatrix->makeSpecialColumnCopy();
    309                 }
    310 #ifdef JJF_ZERO
    311                 if (clpMatrix) {
    312                     int numberRows = clpMatrix->getNumRows();
    313                     int numberColumns = clpMatrix->getNumCols();
    314                     double * elements = clpMatrix->getMutableElements();
    315                     const int * row = clpMatrix->getIndices();
    316                     const CoinBigIndex * columnStart = clpMatrix->getVectorStarts();
    317                     const int * columnLength = clpMatrix->getVectorLengths();
    318                     double * smallest = new double [numberRows];
    319                     double * largest = new double [numberRows];
    320                     char * flag = new char [numberRows];
    321                     CoinZeroN(flag, numberRows);
    322                     for (int i = 0; i < numberRows; i++) {
    323                         smallest[i] = COIN_DBL_MAX;
    324                         largest[i] = 0.0;
    325                     }
    326                     for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
    327                         bool isInteger = modelC->isInteger(iColumn);
    328                         CoinBigIndex j;
    329                         for (j = columnStart[iColumn];
    330                                 j < columnStart[iColumn] + columnLength[iColumn]; j++) {
    331                             int iRow = row[j];
    332                             double value = fabs(elements[j]);
    333                             if (!isInteger)
    334                                 flag[iRow] = 1;
    335                             smallest[iRow] = CoinMin(smallest[iRow], value);
    336                             largest[iRow] = CoinMax(largest[iRow], value);
    337                         }
    338                     }
    339                     double * rowLower = modelC->rowLower();
    340                     double * rowUpper = modelC->rowUpper();
    341                     bool changed = false;
    342                     for (int i = 0; i < numberRows; i++) {
    343                         if (flag[i] && smallest[i] > 10.0 && false) {
    344                             smallest[i] = 1.0 / smallest[i];
    345                             if (rowLower[i] > -1.0e20)
    346                                 rowLower[i] *= smallest[i];
    347                             if (rowUpper[i] < 1.0e20)
    348                                 rowUpper[i] *= smallest[i];
    349                             changed = true;
    350                         } else {
    351                             smallest[i] = 0.0;
    352                         }
    353                     }
    354                     if (changed) {
    355                         printf("SCALED\n");
    356                         for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
    357                             CoinBigIndex j;
    358                             for (j = columnStart[iColumn];
    359                                     j < columnStart[iColumn] + columnLength[iColumn]; j++) {
    360                                 int iRow = row[j];
    361                                 if (smallest[iRow])
    362                                     elements[j] *= smallest[iRow];
    363                             }
    364                         }
    365                     }
    366                     delete [] smallest;
    367                     delete [] largest;
    368                     delete [] flag;
    369                 }
    370 #endif
    371                 model->checkModel();
    372                 modelC->tightenPrimalBounds(0.0, 0, true);
    373                 model->initialSolve();
    374                 if (modelC->dualBound() == 1.0e10) {
    375                     // user did not set - so modify
    376                     // get largest scaled away from bound
    377                     ClpSimplex temp = *modelC;
    378                     temp.dual(0, 7);
    379                     double largestScaled = 1.0e-12;
    380                     double largest = 1.0e-12;
    381                     int numberRows = temp.numberRows();
    382                     const double * rowPrimal = temp.primalRowSolution();
    383                     const double * rowLower = temp.rowLower();
    384                     const double * rowUpper = temp.rowUpper();
    385                     const double * rowScale = temp.rowScale();
    386                     int iRow;
    387                     for (iRow = 0; iRow < numberRows; iRow++) {
    388                         double value = rowPrimal[iRow];
    389                         double above = value - rowLower[iRow];
    390                         double below = rowUpper[iRow] - value;
    391                         if (above < 1.0e12) {
    392                             largest = CoinMax(largest, above);
    393                         }
    394                         if (below < 1.0e12) {
    395                             largest = CoinMax(largest, below);
    396                         }
    397                         if (rowScale) {
    398                             double multiplier = rowScale[iRow];
    399                             above *= multiplier;
    400                             below *= multiplier;
    401                         }
    402                         if (above < 1.0e12) {
    403                             largestScaled = CoinMax(largestScaled, above);
    404                         }
    405                         if (below < 1.0e12) {
    406                             largestScaled = CoinMax(largestScaled, below);
    407                         }
    408                     }
    409 
    410                     int numberColumns = temp.numberColumns();
    411                     const double * columnPrimal = temp.primalColumnSolution();
    412                     const double * columnLower = temp.columnLower();
    413                     const double * columnUpper = temp.columnUpper();
    414                     const double * columnScale = temp.columnScale();
    415                     int iColumn;
    416                     for (iColumn = 0; iColumn < numberColumns; iColumn++) {
    417                         double value = columnPrimal[iColumn];
    418                         double above = value - columnLower[iColumn];
    419                         double below = columnUpper[iColumn] - value;
    420                         if (above < 1.0e12) {
    421                             largest = CoinMax(largest, above);
    422                         }
    423                         if (below < 1.0e12) {
    424                             largest = CoinMax(largest, below);
    425                         }
    426                         if (columnScale) {
    427                             double multiplier = 1.0 / columnScale[iColumn];
    428                             above *= multiplier;
    429                             below *= multiplier;
    430                         }
    431                         if (above < 1.0e12) {
    432                             largestScaled = CoinMax(largestScaled, above);
    433                         }
    434                         if (below < 1.0e12) {
    435                             largestScaled = CoinMax(largestScaled, below);
    436                         }
    437                     }
    438                     std::cout << "Largest (scaled) away from bound " << largestScaled
    439                               << " unscaled " << largest << std::endl;
    440 #ifdef JJF_ZERO
    441                     modelC->setDualBound(CoinMax(1.0001e8,
    442                                                  CoinMin(1000.0*largestScaled, 1.00001e10)));
    443 #else
    444                     modelC->setDualBound(CoinMax(1.0001e9,
    445                                                  CoinMin(1000.0*largestScaled, 1.0001e10)));
    446 #endif
    447                 }
    448             }
    449             model->setMinimumDrop(CoinMin(5.0e-2,
    450                                           fabs(model->getMinimizationObjValue())*1.0e-3 + 1.0e-4));
    451             if (CoinAbs(model->getMaximumCutPassesAtRoot()) <= 100) {
    452                 if (model->getNumCols() < 500) {
    453                     model->setMaximumCutPassesAtRoot(-100); // always do 100 if possible
    454                 } else if (model->getNumCols() < 5000) {
    455                     model->setMaximumCutPassesAtRoot(100); // use minimum drop
    456                 } else {
    457                     model->setMaximumCutPassesAtRoot(20);
    458                 }
    459             }
    460             // If defaults then increase trust for small models
    461             if (model->numberStrong() == 5 && model->numberBeforeTrust() == 10) {
    462                 int numberColumns = model->getNumCols();
    463                 if (numberColumns <= 50) {
    464                     model->setNumberBeforeTrust(1000);
    465                 } else if (numberColumns <= 100) {
    466                     model->setNumberBeforeTrust(100);
    467                 } else if (numberColumns <= 300) {
    468                     model->setNumberBeforeTrust(50);
    469                 }
    470             }
    471             //if (model->getNumCols()>=500) {
    472             // switch off Clp stuff
    473             //model->setFastNodeDepth(-1);
    474             //}
    475             if (model->getNumCols() == -2756) {
    476                 // p2756
    477                 std::string problemName ;
    478                 model->solver()->getStrParam(OsiProbName, problemName) ;
    479                 model->solver()->activateRowCutDebugger(problemName.c_str()) ;
    480             }
    481             if (model->getNumCols() == -201) {
    482                 // p201
    483                 std::string problemName ;
    484                 model->solver()->getStrParam(OsiProbName, problemName) ;
    485                 model->solver()->activateRowCutDebugger(problemName.c_str()) ;
    486             }
    487             if (model->getNumCols() == -104) {
    488                 // bell5
    489                 std::string problemName ;
    490                 model->solver()->getStrParam(OsiProbName, problemName) ;
    491                 model->solver()->activateRowCutDebugger(problemName.c_str()) ;
    492             }
    493             if (model->getNumCols() == -548 && model->getNumRows() == 176) {
    494                 // p0548
    495                 std::string problemName ;
    496                 model->solver()->getStrParam(OsiProbName, problemName) ;
    497                 model->solver()->activateRowCutDebugger(problemName.c_str()) ;
    498             }
    499             if (model->getNumCols() == -160) {
    500                 // misc03
    501                 std::string problemName ;
    502                 model->solver()->getStrParam(OsiProbName, problemName) ;
    503                 model->solver()->activateRowCutDebugger(problemName.c_str()) ;
    504             }
    505             if (model->getNumCols() == -353) {
    506                 // blend2
    507                 std::string problemName ;
    508                 model->solver()->getStrParam(OsiProbName, problemName) ;
    509                 model->solver()->activateRowCutDebugger(problemName.c_str()) ;
    510             }
    511             if (model->getNumCols() == -100 && model->getNumRows() == 21) {
    512                 // enigma
    513                 std::string problemName ;
    514                 model->solver()->getStrParam(OsiProbName, problemName) ;
    515                 model->solver()->activateRowCutDebugger(problemName.c_str()) ;
    516             }
    517             if (model->getNumCols() == -1541) {
    518                 // qnet1
    519                 std::string problemName ;
    520                 model->solver()->getStrParam(OsiProbName, problemName) ;
    521                 model->solver()->activateRowCutDebugger(problemName.c_str()) ;
    522             }
    523             if (model->getNumCols() == -10724) {
    524                 // mitre
    525                 std::string problemName ;
    526                 model->solver()->getStrParam(OsiProbName, problemName) ;
    527                 model->solver()->activateRowCutDebugger(problemName.c_str()) ;
    528             }
    529             if (model->getNumCols() == -1224) {
    530                 //PUSH_MPS("gesa2",1392,1224,25779856.372,25476489.678,7);
    531                 // gesa2
    532                 std::string problemName ;
    533                 model->solver()->getStrParam(OsiProbName, problemName) ;
    534                 model->solver()->activateRowCutDebugger(problemName.c_str()) ;
    535             }
    536             if (model->getNumCols() == -1224 && model->getNumRows() < 1380) {
    537                 //PUSH_MPS("gesa2_o",1248,1224,25779856.372,25476489.678,1);
    538                 // gesa2_o
    539                 std::string problemName ;
    540                 model->solver()->getStrParam(OsiProbName, problemName) ;
    541                 model->solver()->activateRowCutDebugger(problemName.c_str()) ;
    542             }
    543             if (model->getNumCols() == -1152 && model->getNumRows() == 1368) {
    544                 //PUSH_MPS("gesa3",1368,1152,27991042.648,27833632.451,7);
    545                 // gesa3
    546                 std::string problemName ;
    547                 model->solver()->getStrParam(OsiProbName, problemName) ;
    548                 model->solver()->activateRowCutDebugger(problemName.c_str()) ;
    549             }
    550             if (model->getNumCols() == -1152 && model->getNumRows() == 1224) {
    551                 //PUSH_MPS("gesa3_o",1224,1152,27991042.648,27833632.451,7);
    552                 // gesa3
    553                 std::string problemName ;
    554                 model->solver()->getStrParam(OsiProbName, problemName) ;
    555                 model->solver()->activateRowCutDebugger(problemName.c_str()) ;
    556             }
    557             if (model->getNumCols() == -282) {
    558                 //PUSH_MPS("p0282",241,282,258411,176867.50,7);
    559                 // p0282
    560                 std::string problemName ;
    561                 model->solver()->getStrParam(OsiProbName, problemName) ;
    562                 model->solver()->activateRowCutDebugger(problemName.c_str()) ;
    563             }
    564             if (model->getNumCols() == -141) {
    565                 // egout
    566                 std::string problemName ;
    567                 model->solver()->getStrParam(OsiProbName, problemName) ;
    568                 model->solver()->activateRowCutDebugger(problemName.c_str()) ;
    569             }
    570             if (model->getNumCols() == -378) {
    571                 // vpm2
    572                 std::string problemName ;
    573                 model->solver()->getStrParam(OsiProbName, problemName) ;
    574                 model->solver()->activateRowCutDebugger(problemName.c_str()) ;
    575             }
    576             if (model->getNumCols() == -240 && model->getNumRows() == 246) {
    577                 // pp08aCUTS
    578                 std::string problemName ;
    579                 model->solver()->getStrParam(OsiProbName, problemName) ;
    580                 model->solver()->activateRowCutDebugger(problemName.c_str()) ;
    581             }
    582             if (model->getNumCols() == -240 && model->getNumRows() == 136) {
    583                 // pp08a
    584                 std::string problemName ;
    585                 model->solver()->getStrParam(OsiProbName, problemName) ;
    586                 model->solver()->activateRowCutDebugger(problemName.c_str()) ;
    587             }
    588             if (model->getNumCols() == -1372 && model->getNumRows() == 4944) {
    589                 // seymour1
    590                 std::string problemName ;
    591                 model->solver()->getStrParam(OsiProbName, problemName) ;
    592                 model->solver()->activateRowCutDebugger(problemName.c_str()) ;
    593             }
    594             setCutAndHeuristicOptions(*model);
    595             if (si) {
    596 #ifdef CLP_MULTIPLE_FACTORIZATIONS
    597                 if (!modelC->factorization()->isDenseOrSmall()) {
    598                     int denseCode = stuff ? static_cast<int> (stuff[4]) : -1;
    599                     int smallCode = stuff ? static_cast<int> (stuff[10]) : -1;
    600                     if (stuff && stuff[8] >= 1) {
    601                         if (denseCode < 0)
    602                             denseCode = 40;
    603                         if (smallCode < 0)
    604                             smallCode = 40;
    605                     }
    606                     if (denseCode > 0)
    607                         modelC->factorization()->setGoDenseThreshold(denseCode);
    608                     if (smallCode > 0)
    609                         modelC->factorization()->setGoSmallThreshold(smallCode);
    610                     if (denseCode >= modelC->numberRows()) {
    611                         //printf("problem going dense\n");
    612                         //modelC->factorization()->goDenseOrSmall(modelC->numberRows());
    613                     }
    614                 }
    615 #endif
    616                 if (stuff && stuff[8] >= 1) {
    617                     if (modelC->numberColumns() + modelC->numberRows() <= 10000 &&
    618                             model->fastNodeDepth() == -1)
    619                         model->setFastNodeDepth(-9);
    620                 }
    621             }
    622             //OsiObject * obj = new CbcBranchToFixLots(model,0.3,0.0,3,3000003);
    623             //model->addObjects(1,&obj);
    624             //delete obj;
    625             model->branchAndBound();
    626 #ifdef CLP_FACTORIZATION_INSTRUMENT
    627             double facTime = factorization_instrument(0);
    628             printf("Factorization %g seconds\n",
    629                    facTime);
    630             timeTakenFac += facTime;
    631 #endif
    632 
    633             double timeOfSolution = CoinCpuTime() + CoinCpuTimeJustChildren() - startTime;
    634             // Print more statistics
    635             std::cout << "Cuts at root node changed objective from " << model->getContinuousObjective()
    636                       << " to " << model->rootObjectiveAfterCuts() << std::endl;
    637             int numberGenerators = model->numberCutGenerators();
    638             for (int iGenerator = 0; iGenerator < numberGenerators; iGenerator++) {
    639                 CbcCutGenerator * generator = model->cutGenerator(iGenerator);
    640 #ifndef CLP_INVESTIGATE
    641                 CglImplication * implication = dynamic_cast<CglImplication*>(generator->generator());
    642                 if (implication)
    643                     continue;
    644 #endif
    645                 std::cout << generator->cutGeneratorName() << " was tried "
    646                           << generator->numberTimesEntered() << " times and created "
    647                           << generator->numberCutsInTotal() << " cuts of which "
    648                           << generator->numberCutsActive() << " were active after adding rounds of cuts";
    649                 if (generator->timing())
    650                     std::cout << " ( " << generator->timeInCutGenerator() << " seconds)" << std::endl;
    651                 else
    652                     std::cout << std::endl;
    653             }
    654             printf("%d solutions found by heuristics\n",
    655                    model->getNumberHeuristicSolutions());
    656             for (int iHeuristic = 0; iHeuristic < model->numberHeuristics(); iHeuristic++) {
    657                 CbcHeuristic * heuristic = model->heuristic(iHeuristic);
    658                 if (heuristic->numRuns()) {
    659                     // Need to bring others inline
    660                     char generalPrint[1000];
    661                     sprintf(generalPrint, "%s was tried %d times out of %d and created %d solutions\n",
    662                             heuristic->heuristicName(),
    663                             heuristic->numRuns(),
    664                             heuristic->numCouldRun(),
    665                             heuristic->numberSolutionsFound());
    666                     std::cout << generalPrint << std::endl;
    667                 }
    668             }
    669             if (!model->status()) {
    670                 double soln = model->getObjValue();
    671                 double tolerance = CoinMax(1.0e-5, 1.0e-5 * CoinMin(fabs(soln), fabs(objValue[m])));
    672                 //CoinRelFltEq eq(1.0e-3) ;
    673                 if (fabs(soln - objValue[m]) < tolerance) {
    674                     std::cout
    675                         << "cbc_clp" << " "
    676                         << soln << " = " << objValue[m] << " ; okay";
    677                     numProbSolved++;
    678                 } else  {
    679                     std::cout << "cbc_clp" << " " << soln << " != " << objValue[m]
    680                               << "; error=" << fabs(objValue[m] - soln);
    681                     numberFailures++;
    682                     //#ifdef COIN_DEVELOP
    683                     //abort();
    684                     //#endif
    685                 }
    686             } else {
    687                 std::cout << "cbc_clp error; too many nodes" ;
    688             }
    689             timeTaken += timeOfSolution;
    690             std::cout << " - took " << timeOfSolution << " seconds.(" <<
    691                       model->getNodeCount() << " / " << model->getIterationCount() <<
    692                       " ) subtotal " << timeTaken
    693                       << " (" << mpsName[m] << ")" << std::endl;
    694             delete model;
    695         }
    696     }   // end main loop on MPS problem
    697     int returnCode = 0;
     671      std::cout
     672        << "cbc_clp (" << mpsName[m] << ") status not optimal; "
     673        << "assuming too many nodes" ;
     674    }
     675    timeTaken += timeOfSolution;
    698676    std::cout
    699         << "cbc_clp"
    700         << " solved "
    701         << numProbSolved
    702         << " out of "
    703         << numberAttempts;
    704     int numberOnNodes = numberAttempts - numProbSolved - numberFailures;
     677      << " -- (" << model->getNodeCount() << " n / "
     678      << model->getIterationCount() << " i / "
     679      << timeOfSolution << " s) (subtotal " << timeTaken << " s)"
     680      << std::endl;
     681    delete model;
     682  }
     683/*
     684  End main loop on MPS problems. Print a summary and calculate the return
     685  value.
     686*/
     687  int returnCode = 0;
     688  std::cout
     689    << "cbc_clp solved " << numProbSolved << " out of " << numberAttempts;
     690  int numberOnNodes = numberAttempts-numProbSolved-numberFailures;
     691  if (numberFailures || numberOnNodes) {
     692    if (numberOnNodes) {
     693      std::cout << " (" << numberOnNodes << " stopped on nodes)";
     694      returnCode = numberOnNodes;
     695    }
     696    if (numberFailures) {
     697      std::cout << " (" << numberFailures << " gave bad answer!)";
     698      returnCode += 100*numberFailures;
     699    }
     700  }
     701  std::cout
     702    << " and took " << timeTaken << " seconds." << std::endl;
     703
     704  if (testSwitch == -2) {
    705705    if (numberFailures || numberOnNodes) {
    706         if (numberOnNodes) {
    707             std::cout << " (" << numberOnNodes << " stopped on nodes)";
    708             returnCode = numberOnNodes;
    709         }
    710         if (numberFailures) {
    711             std::cout << " (" << numberFailures << " gave bad answer!)";
    712             returnCode += 100 * numberFailures;
    713         }
    714     }
    715     std::cout << " and took "
    716               << timeTaken
    717               << " seconds."
    718               << std::endl;
    719     if (testSwitch == -2) {
    720         if (numberFailures || numberOnNodes) {
    721             printf("****** Unit Test failed\n");
    722             fprintf(stderr, "****** Unit Test failed\n");
    723         } else {
    724             fprintf(stderr, "Unit Test succeeded\n");
    725         }
    726     }
    727 #ifdef CLP_FACTORIZATION_INSTRUMENT
    728     printf("Total factorization time %g seconds\n",
    729            timeTakenFac);
    730 #endif
    731     return returnCode;
     706      std::cout << "****** Unit Test failed." << std::endl ;
     707      std::cerr << "****** Unit Test failed." << std::endl ;
     708    } else {
     709      std::cerr << "****** Unit Test succeeded." << std::endl ;
     710    }
     711  }
     712# ifdef CLP_FACTORIZATION_INSTRUMENT
     713  std::cout
     714    << "Total factorization time " << timeTakenFac << "seconds." << std::endl ;
     715# endif
     716  return (returnCode) ;
    732717}
    733718
Note: See TracChangeset for help on using the changeset viewer.