Changeset 1878


Ignore:
Timestamp:
Aug 30, 2012 11:43:19 AM (7 years ago)
Author:
forrest
Message:

minor changes to implement Aboca

Location:
trunk/Clp/src
Files:
56 added
23 edited

Legend:

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

    r1856 r1878  
    4242#define CBC_THREAD
    4343#endif
     44#endif
     45#if CLP_HAS_ABC
     46#include "AbcCommon.hpp"
    4447#endif
    4548static bool doPrinting = true;
     
    13121315          CbcOrClpParam("-", "From stdin",
    13131316                        CLP_PARAM_ACTION_STDIN, 3, 0);
     1317#ifdef ABC_INHERIT
     1318      parameters[numberParameters++] =
     1319          CbcOrClpParam("abc", "Whether to visit Aboca",
     1320                        "off", CLP_PARAM_STR_ABCWANTED, 7, 0);
     1321     parameters[numberParameters-1].append("one");
     1322     parameters[numberParameters-1].append("two");
     1323     parameters[numberParameters-1].append("three");
     1324     parameters[numberParameters-1].append("four");
     1325     parameters[numberParameters-1].append("five");
     1326     parameters[numberParameters-1].append("six");
     1327     parameters[numberParameters-1].append("seven");
     1328     parameters[numberParameters-1].append("eight");
     1329     parameters[numberParameters-1].append("on");
     1330     parameters[numberParameters-1].append("decide");
     1331     parameters[numberParameters-1].setLonghelp
     1332     (
     1333          "Decide whether to use A Basic Optimization Code (Accelerated?) \
     1334and whether to try going parallel!"
     1335     );
     1336#endif
    13141337     parameters[numberParameters++] =
    13151338          CbcOrClpParam("allC!ommands", "Whether to print less used commands",
     
    15581581     parameters[numberParameters-1].append("so!low_halim");
    15591582     parameters[numberParameters-1].append("lots");
    1560      //  parameters[numberParameters-1].append("4");
    1561      //  parameters[numberParameters-1].append("5");
     1583#ifdef ABC_INHERIT
     1584     parameters[numberParameters-1].append("dual");
     1585     parameters[numberParameters-1].append("dw");
     1586     parameters[numberParameters-1].append("idiot");
     1587#endif
    15621588     parameters[numberParameters-1].setLonghelp
    15631589     (
     
    19852011     parameters[numberParameters-1].setLonghelp
    19862012     (
     2013#ifndef ABC_INHERIT
    19872014          "The default is to use the normal CoinFactorization, but \
    19882015other choices are a dense one, osl's or one designed for small problems."
     2016#else
     2017          "Normally the default is to use the normal CoinFactorization, but \
     2018other choices are a dense one, osl's or one designed for small problems. \
     2019However if at Aboca then the default is CoinAbcFactorization and other choices are \
     2020a dense one, one designed for small problems or if enabled a long factorization."
     2021#endif
    19892022     );
    19902023     parameters[numberParameters++] =
     
    28212854     parameters[numberParameters-1].append("rhs!ranging");
    28222855     parameters[numberParameters-1].append("objective!ranging");
     2856     parameters[numberParameters-1].append("stats");
    28232857     parameters[numberParameters-1].setLonghelp
    28242858     (
     
    31713205     );
    31723206#else
    3173      // allow solve as synonym for dual
    3174      parameters[numberParameters++] =
    3175           CbcOrClpParam("solv!e", "Solve problem using dual simplex",
    3176                         CBC_PARAM_ACTION_BAB);
     3207     // allow solve as synonym for possible dual
     3208     parameters[numberParameters++] =
     3209       CbcOrClpParam("solv!e", "Solve problem using dual simplex (probably)",
     3210                        CLP_PARAM_ACTION_EITHERSIMPLEX);
    31773211     parameters[numberParameters-1].setLonghelp
    31783212     (
     
    33223356          "For making equality cliques this is minimumsize.  Also for adding \
    33233357integer slacks.  May be used for more later \
    3324 If <1000 that is what it does.  If <1000000 - numberPasses is (value/1000)-1 and tune is tune %1000. \
     3358If <10000 that is what it does.  If <1000000 - numberPasses is (value/10000)-1 and tune is tune %10000. \
    33253359If >= 1000000! - numberPasses is (value/1000000)-1 and tune is tune %1000000.  In this case if tune is now still >=10000 \
    33263360numberPassesPerInnerLoop is changed from 10 to (tune-10000)-1 and tune becomes tune % 10000!!!!! - happy? - \
  • trunk/Clp/src/CbcOrClpParam.hpp

    r1852 r1878  
    160160     CLP_PARAM_STR_ALLCOMMANDS,
    161161     CLP_PARAM_STR_TIME_MODE,
     162     CLP_PARAM_STR_ABCWANTED,
    162163
    163164     CBC_PARAM_STR_NODESTRATEGY = 251,
  • trunk/Clp/src/ClpCholeskyBase.cpp

    r1732 r1878  
    33003300     bool endClique = false;
    33013301     int lastRow = 0;
    3302      CoinBigIndex cliquePointer = 0;
    33033302     int nextRow2 = -1;
    33043303
     
    33963395               // initialize new clique
    33973396               lastRow = iRow;
    3398                cliquePointer = choleskyStart_[iRow];
    33993397          }
    34003398          // for each column L[*,kRow] that affects L[*,iRow]
  • trunk/Clp/src/ClpCholeskyDense.cpp

    r1726 r1878  
    1414#include "ClpMessage.hpp"
    1515#include "ClpQuadraticObjective.hpp"
     16#if CLP_HAS_ABC
     17#include "CoinAbcCommon.hpp"
     18#endif
     19#if CLP_HAS_ABC<3
     20#undef cilk_for
     21#undef cilk_spawn
     22#undef cilk_sync
     23#define cilk_for for
     24#define cilk_spawn
     25#define cilk_sync
     26#endif
    1627
    1728/*#############################################################################*/
     
    695706          int nb = number_blocks((nLeft + 1) >> 1);
    696707          int nLeft2 = number_rows(nb);
    697           ClpCholeskyCtriRec(thisStruct, aTri, nThis, aUnder, diagonal, work, nLeft2, iBlock, jBlock, numberBlocks);
     708          cilk_spawn ClpCholeskyCtriRec(thisStruct, aTri, nThis, aUnder, diagonal, work, nLeft2, iBlock, jBlock, numberBlocks);
    698709          ClpCholeskyCtriRec(thisStruct, aTri, nThis, aUnder + number_entries(nb), diagonal, work, nLeft - nLeft2,
    699710                             iBlock + nb, jBlock, numberBlocks);
     711          cilk_sync;
    700712     } else {
    701713          int nb = number_blocks((nThis + 1) >> 1);
     
    743755          longDouble * aother;
    744756          int i;
    745           ClpCholeskyCrecTri(thisStruct, aUnder, nTri2, nDo, iBlock, jBlock, aTri, diagonal, work, numberBlocks);
     757          cilk_spawn ClpCholeskyCrecTri(thisStruct, aUnder, nTri2, nDo, iBlock, jBlock, aTri, diagonal, work, numberBlocks);
    746758          /* and rectangular update */
    747759          i = ((numberBlocks - iBlock) * (numberBlocks - iBlock + 1) -
     
    752764          ClpCholeskyCrecTri(thisStruct, aUnder + number_entries(nb), nTri - nTri2, nDo, iBlock + nb, jBlock,
    753765                             aTri + number_entries(i), diagonal, work, numberBlocks);
     766          cilk_sync;
    754767     }
    755768}
     
    771784          int nb = number_blocks((nUnderK + 1) >> 1);
    772785          int nUnder2 = number_rows(nb);
    773           ClpCholeskyCrecRec(thisStruct, above, nUnder, nUnder2, nDo, aUnder, aOther, work,
     786          cilk_spawn ClpCholeskyCrecRec(thisStruct, above, nUnder, nUnder2, nDo, aUnder, aOther, work,
    774787                             iBlock, jBlock, numberBlocks);
    775788          ClpCholeskyCrecRec(thisStruct, above, nUnder, nUnderK - nUnder2, nDo, aUnder + number_entries(nb),
    776789                             aOther + number_entries(nb), work, iBlock, jBlock, numberBlocks);
     790          cilk_sync;
    777791     } else if (nUnderK <= nDo && nUnder <= nDo) {
    778792          int nb = number_blocks((nDo + 1) >> 1);
     
    791805          int nUnder2 = number_rows(nb);
    792806          int i;
    793           ClpCholeskyCrecRec(thisStruct, above, nUnder2, nUnderK, nDo, aUnder, aOther, work,
     807          cilk_spawn ClpCholeskyCrecRec(thisStruct, above, nUnder2, nUnderK, nDo, aUnder, aOther, work,
    794808                             iBlock, jBlock, numberBlocks);
    795809          i = ((numberBlocks - iBlock) * (numberBlocks - iBlock - 1) -
     
    797811          ClpCholeskyCrecRec(thisStruct, above + number_entries(nb), nUnder - nUnder2, nUnderK, nDo, aUnder,
    798812                             aOther + number_entries(i), work, iBlock + nb, jBlock, numberBlocks);
     813          cilk_sync;
    799814     }
    800815}
  • trunk/Clp/src/ClpMain.cpp

    r1726 r1878  
    3737#endif
    3838
     39#if CLP_HAS_ABC
     40#include "AbcCommon.hpp"
     41#endif
    3942#include "ClpFactorization.hpp"
    4043#include "CoinTime.hpp"
     
    5356#include "CbcOrClpParam.hpp"
    5457#include "CoinSignal.hpp"
     58#ifdef ABC_INHERIT
     59#include "AbcSimplex.hpp"
     60#include "AbcSimplexFactorization.hpp"
     61#include "AbcDualRowSteepest.hpp"
     62#include "AbcDualRowDantzig.hpp"
     63#endif
    5564#ifdef DMALLOC
    5665#include "dmalloc.h"
     
    6372static bool maskMatches(const int * starts, char ** masks,
    6473                        std::string & check);
     74#ifndef ABC_INHERIT
    6575static ClpSimplex * currentModel = NULL;
     76#else
     77static AbcSimplex * currentModel = NULL;
     78#endif
    6679
    6780extern "C" {
     
    8497#endif
    8598
     99#ifndef ABC_INHERIT
    86100int mainTest (int argc, const char *argv[], int algorithm,
    87101              ClpSimplex empty, ClpSolve solveOptions, int switchOff, bool doVector);
     102#else
     103int mainTest (int argc, const char *argv[], int algorithm,
     104              AbcSimplex empty, ClpSolve solveOptions, int switchOff, bool doVector);
     105#endif
    88106static void statistics(ClpSimplex * originalModel, ClpSimplex * model);
    89107static void generateCode(const char * fileName, int type);
     
    127145}
    128146#endif
    129 
     147//#define CILK_TEST
     148#ifdef CILK_TEST
     149static void cilkTest();
     150#endif
    130151int
    131152#if defined(_MSC_VER)
     
    134155main (int argc, const char *argv[])
    135156{
     157#ifdef CILK_TEST
     158  cilkTest();
     159#endif
    136160     // next {} is just to make sure all memory should be freed - for debug
    137161     {
     
    139163          // Set up all non-standard stuff
    140164          //int numberModels=1;
     165#ifndef ABC_INHERIT
    141166          ClpSimplex * models = new ClpSimplex[1];
     167#else
     168          AbcSimplex * models = new AbcSimplex[1];
     169#endif
    142170
    143171          // default action on import
     
    155183          int doSprint = -1;
    156184          // set reasonable defaults
    157           int preSolve = 5;
     185#ifdef ABC_INHERIT
     186#define DEFAULT_PRESOLVE_PASSES 20
     187#else
     188#define DEFAULT_PRESOLVE_PASSES 5
     189#endif
     190          int preSolve = DEFAULT_PRESOLVE_PASSES;
    158191          bool preSolveFile = false;
    159192          models->setPerturbation(50);
     
    248281          //ClpDualRowSteepest steep;
    249282          //models[0].setDualRowPivotAlgorithm(steep);
    250           //models[0].setPrimalTolerance(1.0e-7);
     283#ifdef ABC_INHERIT
     284          models[0].setDualTolerance(1.0e-6);
     285          models[0].setPrimalTolerance(1.0e-6);
     286#endif
    251287          //ClpPrimalColumnSteepest steepP;
    252288          //models[0].setPrimalColumnPivotAlgorithm(steepP);
     
    312348                    }
    313349               }
     350               ClpSimplex * thisModel=models+iModel;
    314351               if (iParam < numberParameters && !numberQuery) {
    315352                    // found
     
    431468                         double value = CoinReadGetDoubleField(argc, argv, &valid);
    432469                         if (!valid) {
    433                               parameters[iParam].setDoubleParameter(models + iModel, value);
     470                              parameters[iParam].setDoubleParameter(thisModel, value);
    434471                         } else if (valid == 1) {
    435472                              std::cout << " is illegal for double parameter " << parameters[iParam].name() << " value remains " <<
     
    465502                              else if (parameters[iParam].type() == CLP_PARAM_INT_VERBOSE)
    466503                                   verbose = value;
    467                               parameters[iParam].setIntParameter(models + iModel, value);
     504                              parameters[iParam].setIntParameter(thisModel, value);
    468505                         } else if (valid == 1) {
    469506                              std::cout << " is illegal for integer parameter " << parameters[iParam].name() << " value remains " <<
     
    491528                              switch (type) {
    492529                              case CLP_PARAM_STR_DIRECTION:
    493                                    if (action == 0)
     530                                if (action == 0) {
    494531                                        models[iModel].setOptimizationDirection(1);
    495                                    else if (action == 1)
     532#ifdef ABC_INHERIT
     533                                        thisModel->setOptimizationDirection(1);
     534#endif
     535                                }  else if (action == 1) {
    496536                                        models[iModel].setOptimizationDirection(-1);
    497                                    else
     537#ifdef ABC_INHERIT
     538                                        thisModel->setOptimizationDirection(-1);
     539#endif
     540                                }  else {
    498541                                        models[iModel].setOptimizationDirection(0);
     542#ifdef ABC_INHERIT
     543                                        thisModel->setOptimizationDirection(0);
     544#endif
     545                                }
    499546                                   break;
    500547                              case CLP_PARAM_STR_DUALPIVOT:
    501548                                   if (action == 0) {
    502549                                        ClpDualRowSteepest steep(3);
    503                                         models[iModel].setDualRowPivotAlgorithm(steep);
     550                                        thisModel->setDualRowPivotAlgorithm(steep);
     551#ifdef ABC_INHERIT
     552                                        AbcDualRowSteepest steep2(3);
     553                                        models[iModel].setDualRowPivotAlgorithm(steep2);
     554#endif
    504555                                   } else if (action == 1) {
    505556                                        //ClpDualRowDantzig dantzig;
    506                                         ClpDualRowSteepest dantzig(5);
    507                                         models[iModel].setDualRowPivotAlgorithm(dantzig);
     557                                        ClpDualRowDantzig dantzig;
     558                                        thisModel->setDualRowPivotAlgorithm(dantzig);
     559#ifdef ABC_INHERIT
     560                                        AbcDualRowDantzig dantzig2;
     561                                        models[iModel].setDualRowPivotAlgorithm(dantzig2);
     562#endif
    508563                                   } else if (action == 2) {
    509564                                        // partial steep
    510565                                        ClpDualRowSteepest steep(2);
    511                                         models[iModel].setDualRowPivotAlgorithm(steep);
     566                                        thisModel->setDualRowPivotAlgorithm(steep);
     567#ifdef ABC_INHERIT
     568                                        AbcDualRowSteepest steep2(2);
     569                                        models[iModel].setDualRowPivotAlgorithm(steep2);
     570#endif
    512571                                   } else {
    513572                                        ClpDualRowSteepest steep;
    514                                         models[iModel].setDualRowPivotAlgorithm(steep);
     573                                        thisModel->setDualRowPivotAlgorithm(steep);
     574#ifdef ABC_INHERIT
     575                                        AbcDualRowSteepest steep2;
     576                                        models[iModel].setDualRowPivotAlgorithm(steep2);
     577#endif
    515578                                   }
    516579                                   break;
     
    518581                                   if (action == 0) {
    519582                                        ClpPrimalColumnSteepest steep(3);
    520                                         models[iModel].setPrimalColumnPivotAlgorithm(steep);
     583                                        thisModel->setPrimalColumnPivotAlgorithm(steep);
    521584                                   } else if (action == 1) {
    522585                                        ClpPrimalColumnSteepest steep(0);
    523                                         models[iModel].setPrimalColumnPivotAlgorithm(steep);
     586                                        thisModel->setPrimalColumnPivotAlgorithm(steep);
    524587                                   } else if (action == 2) {
    525588                                        ClpPrimalColumnDantzig dantzig;
    526                                         models[iModel].setPrimalColumnPivotAlgorithm(dantzig);
     589                                        thisModel->setPrimalColumnPivotAlgorithm(dantzig);
    527590                                   } else if (action == 3) {
    528591                                        ClpPrimalColumnSteepest steep(4);
    529                                         models[iModel].setPrimalColumnPivotAlgorithm(steep);
     592                                        thisModel->setPrimalColumnPivotAlgorithm(steep);
    530593                                   } else if (action == 4) {
    531594                                        ClpPrimalColumnSteepest steep(1);
    532                                         models[iModel].setPrimalColumnPivotAlgorithm(steep);
     595                                        thisModel->setPrimalColumnPivotAlgorithm(steep);
    533596                                   } else if (action == 5) {
    534597                                        ClpPrimalColumnSteepest steep(2);
    535                                         models[iModel].setPrimalColumnPivotAlgorithm(steep);
     598                                        thisModel->setPrimalColumnPivotAlgorithm(steep);
    536599                                   } else if (action == 6) {
    537600                                        ClpPrimalColumnSteepest steep(10);
    538                                         models[iModel].setPrimalColumnPivotAlgorithm(steep);
     601                                        thisModel->setPrimalColumnPivotAlgorithm(steep);
    539602                                   }
    540603                                   break;
    541604                              case CLP_PARAM_STR_SCALING:
    542                                    models[iModel].scaling(action);
     605                                   thisModel->scaling(action);
    543606                                   break;
    544607                              case CLP_PARAM_STR_AUTOSCALE:
    545                                    models[iModel].setAutomaticScaling(action != 0);
     608                                   thisModel->setAutomaticScaling(action != 0);
    546609                                   break;
    547610                              case CLP_PARAM_STR_SPARSEFACTOR:
    548                                    models[iModel].setSparseFactorization((1 - action) != 0);
     611                                   thisModel->setSparseFactorization((1 - action) != 0);
    549612                                   break;
    550613                              case CLP_PARAM_STR_BIASLU:
    551                                    models[iModel].factorization()->setBiasLU(action);
     614                                   thisModel->factorization()->setBiasLU(action);
    552615                                   break;
    553616                              case CLP_PARAM_STR_PERTURBATION:
    554617                                   if (action == 0)
    555                                         models[iModel].setPerturbation(50);
     618                                        thisModel->setPerturbation(50);
    556619                                   else
    557                                         models[iModel].setPerturbation(100);
     620                                        thisModel->setPerturbation(100);
    558621                                   break;
    559622                              case CLP_PARAM_STR_ERRORSALLOWED:
    560623                                   allowImportErrors = action;
     624                                   break;
     625                              case CLP_PARAM_STR_ABCWANTED:
     626                                   models[iModel].setAbcState(action);
    561627                                   break;
    562628                              case CLP_PARAM_STR_INTPRINT:
     
    568634                              case CLP_PARAM_STR_PRESOLVE:
    569635                                   if (action == 0)
    570                                         preSolve = 5;
     636                                        preSolve = DEFAULT_PRESOLVE_PASSES;
    571637                                   else if (action == 1)
    572638                                        preSolve = 0;
     
    577643                                   break;
    578644                              case CLP_PARAM_STR_PFI:
    579                                    models[iModel].factorization()->setForrestTomlin(action == 0);
     645                                   thisModel->factorization()->setForrestTomlin(action == 0);
    580646                                   break;
    581647                              case CLP_PARAM_STR_FACTORIZATION:
    582648                                   models[iModel].factorization()->forceOtherFactorization(action);
     649#ifdef ABC_INHERIT
     650                                   thisModel->factorization()->forceOtherFactorization(action);
     651#endif
    583652                                   break;
    584653                              case CLP_PARAM_STR_CRASH:
     
    590659                              case CLP_PARAM_STR_MESSAGES:
    591660                                   models[iModel].messageHandler()->setPrefix(action != 0);
     661#ifdef ABC_INHERIT
     662                                   thisModel->messageHandler()->setPrefix(action != 0);
     663#endif
    592664                                   break;
    593665                              case CLP_PARAM_STR_CHOLESKY:
     
    623695                         case CBC_PARAM_ACTION_BAB:
    624696                              if (goodModels[iModel]) {
     697                                if (type==CLP_PARAM_ACTION_EITHERSIMPLEX||
     698                                    type==CBC_PARAM_ACTION_BAB)
     699                                  models[iModel].setMoreSpecialOptions(16384|
     700                                                                       models[iModel].moreSpecialOptions());
    625701                                   double objScale =
    626702                                        parameters[whichParam(CLP_PARAM_DBL_OBJSCALE2, numberParameters, parameters)].doubleValue();
     
    648724                                   ClpSolve::SolveType method;
    649725                                   ClpSolve::PresolveType presolveType;
     726                                   ClpSolve solveOptions;
     727#ifndef ABC_INHERIT
    650728                                   ClpSimplex * model2 = models + iModel;
    651                                    ClpSolve solveOptions;
     729#else
     730                                   AbcSimplex * model2 = models + iModel;
     731                                   if (type==CLP_PARAM_ACTION_EITHERSIMPLEX||
     732                                       type==CBC_PARAM_ACTION_BAB)
     733                                     solveOptions.setSpecialOption(3,0); // allow +-1
     734#endif
    652735                                   if (dualize==4) {
    653736                                     solveOptions.setSpecialOption(4, 77);
     
    662745                                             int numberColumns = model2->numberColumns();
    663746                                             int numberRows = model2->numberRows();
     747#ifndef ABC_INHERIT
    664748                                             if (numberRows < 50000 || 5 * numberColumns > numberRows) {
     749#else
     750                                             if (numberRows < 500 || 4 * numberColumns > numberRows) {
     751#endif
    665752                                                  tryIt = false;
    666753                                             } else {
    667754                                                  fractionColumn = 0.1;
    668                                                   fractionRow = 0.1;
     755                                                  fractionRow = 0.3;
    669756                                             }
    670757                                        }
    671758                                        if (tryIt) {
    672                                              model2 = static_cast<ClpSimplexOther *> (model2)->dualOfModel(fractionRow, fractionColumn);
    673                                              if (model2) {
     759                                          ClpSimplex * thisModel=model2;
     760                                             thisModel = static_cast<ClpSimplexOther *> (thisModel)->dualOfModel(fractionRow, fractionColumn);
     761                                             if (thisModel) {
    674762                                                  printf("Dual of model has %d rows and %d columns\n",
    675                                                          model2->numberRows(), model2->numberColumns());
    676                                                   model2->setOptimizationDirection(1.0);
     763                                                         thisModel->numberRows(), thisModel->numberColumns());
     764                                                  thisModel->setOptimizationDirection(1.0);
     765#ifndef ABC_INHERIT
     766                                                  model2=thisModel;
     767#else
     768                                                  int abcState=model2->abcState();
     769                                                  model2=new AbcSimplex(*thisModel);
     770                                                  model2->setAbcState(abcState);
     771                                                  delete thisModel;
     772#endif
    677773                                             } else {
    678                                                   model2 = models + iModel;
     774                                                  thisModel = models + iModel;
    679775                                                  dualize = 0;
    680776                                             }
     
    687783                                   solveOptions.setPresolveActions(presolveOptions);
    688784                                   solveOptions.setSubstitution(substitution);
    689                                    if (preSolve != 5 && preSolve) {
     785                                   if (preSolve != DEFAULT_PRESOLVE_PASSES && preSolve) {
    690786                                        presolveType = ClpSolve::presolveNumber;
    691787                                        if (preSolve < 0) {
     
    718814                                   } else {
    719815                                        method = ClpSolve::useBarrier;
     816#ifdef ABC_INHERIT
     817                                        if (doIdiot > 0)
     818                                             solveOptions.setSpecialOption(1, 2, doIdiot); // dense threshold
     819#endif
    720820                                        if (crossover == 1) {
    721821                                             method = ClpSolve::useBarrierNoCross;
     
    808908#ifdef CLP_MULTIPLE_FACTORIZATIONS
    809909                                   int denseCode = parameters[whichParam(CBC_PARAM_INT_DENSE, numberParameters, parameters)].intValue();
    810                                    model2->factorization()->setGoDenseThreshold(denseCode);
     910                                   if (denseCode!=-1)
     911                                     model2->factorization()->setGoDenseThreshold(denseCode);
    811912                                   int smallCode = parameters[whichParam(CBC_PARAM_INT_SMALLFACT, numberParameters, parameters)].intValue();
    812                                    model2->factorization()->setGoSmallThreshold(smallCode);
     913                                   if (smallCode!=-1)
     914                                     model2->factorization()->setGoSmallThreshold(smallCode);
    813915                                   model2->factorization()->goDenseOrSmall(model2->numberRows());
    814916#endif
     
    820922                                   }
    821923                                   if (dualize) {
    822                                         int returnCode = static_cast<ClpSimplexOther *> (models + iModel)->restoreFromDual(model2);
     924                                     ClpSimplex * thisModel=models+iModel;
     925                                        int returnCode = static_cast<ClpSimplexOther *> (thisModel)->restoreFromDual(model2);
    823926                                        if (model2->status() == 3)
    824927                                             returnCode = 0;
     
    828931                                             // register signal handler
    829932                                             signal(SIGINT, signal_handler);
    830                                              models[iModel].primal(1);
     933                                             thisModel->primal(1);
    831934                                             currentModel = NULL;
    832935                                        }
     936                                        // switch off (user can switch back on)
     937                                        parameters[whichParam(CLP_PARAM_INT_DUALIZE,
     938                                                              numberParameters, parameters)].setIntValue(dualize);
    833939                                   }
    834940                                   if (status >= 0)
     
    10381144                                        goodModels[iModel] = true;
    10391145                                        // sets to all slack (not necessary?)
    1040                                         models[iModel].createStatus();
     1146                                        thisModel->createStatus();
    10411147                                        time2 = CoinCpuTime();
    10421148                                        totalTime += time2 - time1;
     
    12841390                                   }
    12851391                                   if (canOpen) {
    1286                                         int values = models[iModel].readBasis(fileName.c_str());
     1392                                        int values = thisModel->readBasis(fileName.c_str());
    12871393                                        if (values == 0)
    12881394                                             basisHasValues = -1;
     
    15181624                         case CLP_PARAM_ACTION_MAXIMIZE:
    15191625                              models[iModel].setOptimizationDirection(-1);
     1626#ifdef ABC_INHERIT
     1627                              thisModel->setOptimizationDirection(-1);
     1628#endif
    15201629                              break;
    15211630                         case CLP_PARAM_ACTION_MINIMIZE:
    15221631                              models[iModel].setOptimizationDirection(1);
     1632#ifdef ABC_INHERIT
     1633                              thisModel->setOptimizationDirection(1);
     1634#endif
    15231635                              break;
    15241636                         case CLP_PARAM_ACTION_ALLSLACK:
    1525                               models[iModel].allSlackBasis(true);
     1637                              thisModel->allSlackBasis(true);
     1638#ifdef ABC_INHERIT
     1639                              models[iModel].allSlackBasis();
     1640#endif
    15261641                              break;
    15271642                         case CLP_PARAM_ACTION_REVERSE:
     
    16321747                                   std::cerr << "Doing netlib with dual algorithm" << std::endl;
    16331748                                   algorithm = 0;
     1749                                   models[iModel].setMoreSpecialOptions(models[iModel].moreSpecialOptions()|32768);
    16341750                              } else if (type == CLP_PARAM_ACTION_NETLIB_BARRIER) {
    16351751                                   std::cerr << "Doing netlib with barrier algorithm" << std::endl;
     
    16441760                                   // algorithm=6;
    16451761                              } else {
    1646                                    std::cerr << "Doing netlib with primal agorithm" << std::endl;
     1762                                   std::cerr << "Doing netlib with primal algorithm" << std::endl;
    16471763                                   algorithm = 1;
    16481764                              }
    1649                               int specialOptions = models[iModel].specialOptions();
     1765                              //int specialOptions = models[iModel].specialOptions();
    16501766                              models[iModel].setSpecialOptions(0);
    16511767                              ClpSolve solveOptions;
     
    16761792                                   }
    16771793                              }
    1678                               mainTest(nFields, fields, algorithm, models[iModel],
     1794#if FACTORIZATION_STATISTICS
     1795                              {
     1796                                extern int loSizeX;
     1797                                extern int hiSizeX;
     1798                                for (int i=0;i<argc;i++) {
     1799                                  if (!strcmp(argv[i],"-losize")) {
     1800                                    int size=atoi(argv[i+1]);
     1801                                    if (size>0)
     1802                                      loSizeX=size;
     1803                                  }
     1804                                  if (!strcmp(argv[i],"-hisize")) {
     1805                                    int size=atoi(argv[i+1]);
     1806                                    if (size>loSizeX)
     1807                                      hiSizeX=size;
     1808                                  }
     1809                                }
     1810                                if (loSizeX!=-1||hiSizeX!=1000000)
     1811                                  printf("Solving problems %d<= and <%d\n",loSizeX,hiSizeX);
     1812                              }
     1813#endif
     1814                              // for moment then back to models[iModel]
     1815#ifndef ABC_INHERIT
     1816                              int specialOptions = models[iModel].specialOptions();
     1817                              mainTest(nFields, fields, algorithm, *thisModel,
    16791818                                       solveOptions, specialOptions, doVector != 0);
     1819#else
     1820                              //if (!processTune) {
     1821                              //specialOptions=0;
     1822                              //models->setSpecialOptions(models->specialOptions()&~65536);
     1823                              // }
     1824                              mainTest(nFields, fields, algorithm, *models,
     1825                                       solveOptions, 0, doVector != 0);
     1826#endif
    16801827                         }
    16811828                         break;
     
    17001847                                   presolveType = ClpSolve::presolveOff;
    17011848                              solveOptions.setPresolveType(presolveType, 5);
    1702                               mainTest(nFields, fields, algorithm, models[iModel],
     1849#ifndef ABC_INHERIT
     1850                              mainTest(nFields, fields, algorithm, *thisModel,
    17031851                                       solveOptions, specialOptions, doVector != 0);
     1852#else
     1853                              mainTest(nFields, fields, algorithm, *models,
     1854                                       solveOptions, specialOptions, doVector != 0);
     1855#endif
    17041856                         }
    17051857                         break;
     
    19192071                                        }
    19202072                                        fprintf(fp, "Objective value %15.8g\n", objValue);
     2073                                        if (printMode==9) {
     2074                                          // just statistics
     2075                                          int numberRows = models[iModel].numberRows();
     2076                                          double * dualRowSolution = models[iModel].dualRowSolution();
     2077                                          double * primalRowSolution =
     2078                                            models[iModel].primalRowSolution();
     2079                                          double * rowLower = models[iModel].rowLower();
     2080                                          double * rowUpper = models[iModel].rowUpper();
     2081                                          double highestPrimal;
     2082                                          double lowestPrimal;
     2083                                          double highestDual;
     2084                                          double lowestDual;
     2085                                          double largestAway;
     2086                                          int numberAtLower;
     2087                                          int numberAtUpper;
     2088                                          int numberBetween;
     2089                                          highestPrimal=-COIN_DBL_MAX;
     2090                                          lowestPrimal=COIN_DBL_MAX;
     2091                                          highestDual=-COIN_DBL_MAX;
     2092                                          lowestDual=COIN_DBL_MAX;
     2093                                          largestAway=0.0;;
     2094                                          numberAtLower=0;
     2095                                          numberAtUpper=0;
     2096                                          numberBetween=0;
     2097                                          for (int iRow=0;iRow<numberRows;iRow++) {
     2098                                            double primal=primalRowSolution[iRow];
     2099                                            double lower=rowLower[iRow];
     2100                                            double upper=rowUpper[iRow];
     2101                                            double dual=dualRowSolution[iRow];
     2102                                            highestPrimal=CoinMax(highestPrimal,primal);
     2103                                            lowestPrimal=CoinMin(lowestPrimal,primal);
     2104                                            highestDual=CoinMax(highestDual,dual);
     2105                                            lowestDual=CoinMin(lowestDual,dual);
     2106                                            if (primal<lower+1.0e-6) {
     2107                                              numberAtLower++;
     2108                                            } else if (primal>upper-1.0e-6) {
     2109                                              numberAtUpper++;
     2110                                            } else {
     2111                                              numberBetween++;
     2112                                              largestAway=CoinMax(largestAway,
     2113                                                                  CoinMin(primal-lower,upper-primal));
     2114                                            }
     2115                                          }
     2116                                          printf("For rows %d at lower, %d between, %d at upper - lowest %g, highest %g most away %g - highest dual %g lowest %g\n",
     2117                                                 numberAtLower,numberBetween,numberAtUpper,
     2118                                                 lowestPrimal,highestPrimal,largestAway,
     2119                                                 lowestDual,highestDual);
     2120                                          int numberColumns = models[iModel].numberColumns();
     2121                                          double * dualColumnSolution = models[iModel].dualColumnSolution();
     2122                                          double * primalColumnSolution =
     2123                                            models[iModel].primalColumnSolution();
     2124                                          double * columnLower = models[iModel].columnLower();
     2125                                          double * columnUpper = models[iModel].columnUpper();
     2126                                          highestPrimal=-COIN_DBL_MAX;
     2127                                          lowestPrimal=COIN_DBL_MAX;
     2128                                          highestDual=-COIN_DBL_MAX;
     2129                                          lowestDual=COIN_DBL_MAX;
     2130                                          largestAway=0.0;;
     2131                                          numberAtLower=0;
     2132                                          numberAtUpper=0;
     2133                                          numberBetween=0;
     2134                                          for (int iColumn=0;iColumn<numberColumns;iColumn++) {
     2135                                            double primal=primalColumnSolution[iColumn];
     2136                                            double lower=columnLower[iColumn];
     2137                                            double upper=columnUpper[iColumn];
     2138                                            double dual=dualColumnSolution[iColumn];
     2139                                            highestPrimal=CoinMax(highestPrimal,primal);
     2140                                            lowestPrimal=CoinMin(lowestPrimal,primal);
     2141                                            highestDual=CoinMax(highestDual,dual);
     2142                                            lowestDual=CoinMin(lowestDual,dual);
     2143                                            if (primal<lower+1.0e-6) {
     2144                                              numberAtLower++;
     2145                                            } else if (primal>upper-1.0e-6) {
     2146                                              numberAtUpper++;
     2147                                            } else {
     2148                                              numberBetween++;
     2149                                              largestAway=CoinMax(largestAway,
     2150                                                                  CoinMin(primal-lower,upper-primal));
     2151                                            }
     2152                                          }
     2153                                          printf("For columns %d at lower, %d between, %d at upper - lowest %g, highest %g most away %g - highest dual %g lowest %g\n",
     2154                                                 numberAtLower,numberBetween,numberAtUpper,
     2155                                                 lowestPrimal,highestPrimal,largestAway,
     2156                                                 lowestDual,highestDual);
     2157                                          break;
     2158                                        }
    19212159                                        // make fancy later on
    19222160                                        int iRow;
     
    27072945     number = new int[numberColumns+1];
    27082946     memset(number, 0, (numberColumns + 1)*sizeof(int));
     2947     if (model->logLevel() > 3) {
     2948       // get column copy
     2949       CoinPackedMatrix columnCopy = *matrix;
     2950       const int * columnLength = columnCopy.getVectorLengths();
     2951       const int * row = columnCopy.getIndices();
     2952       const CoinBigIndex * columnStart = columnCopy.getVectorStarts();
     2953       const double * element = columnCopy.getElements();
     2954       const double * elementByRow = rowCopy.getElements();
     2955       const int * rowStart = rowCopy.getVectorStarts();
     2956       const int * column = rowCopy.getIndices();
     2957       int nPossibleZeroCost=0;
     2958       int nPossibleNonzeroCost=0;
     2959       for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
     2960         int length = columnLength[iColumn];
     2961         if (columnLower[iColumn]<-1.0e30&&columnUpper[iColumn]>1.0e30) {
     2962           if (length==1) {
     2963             printf("Singleton free %d - cost %g\n",iColumn,objective[iColumn]);
     2964           } else if (length==2) {
     2965             int iRow0=row[columnStart[iColumn]];
     2966             int iRow1=row[columnStart[iColumn]+1];
     2967             double element0=element[columnStart[iColumn]];
     2968             double element1=element[columnStart[iColumn]+1];
     2969             int n0=rowLength[iRow0];
     2970             int n1=rowLength[iRow1];
     2971             printf("Doubleton free %d - cost %g - %g in %srow with %d entries and %g in %srow with %d entries\n",
     2972                    iColumn,objective[iColumn],element0,(rowLower[iRow0]==rowUpper[iRow0]) ? "==" : "",n0,
     2973                    element1,(rowLower[iRow1]==rowUpper[iRow1]) ? "==" : "",n1);
     2974
     2975           }
     2976         }
     2977         if (length==1) {
     2978           int iRow=row[columnStart[iColumn]];
     2979           double value=COIN_DBL_MAX;
     2980           for (int i=rowStart[iRow];i<rowStart[iRow]+rowLength[iRow];i++) {
     2981             int jColumn=column[i];
     2982             if (jColumn!=iColumn) {
     2983               if (value!=elementByRow[i]) {
     2984                 if (value==COIN_DBL_MAX) {
     2985                   value=elementByRow[i];
     2986                 } else {
     2987                   value = -COIN_DBL_MAX;
     2988                   break;
     2989                 }
     2990               }
     2991             }
     2992           }
     2993           if (!objective[iColumn]) {
     2994             if (model->logLevel() > 4)
     2995             printf("Singleton %d with no objective in row with %d elements - rhs %g,%g\n",iColumn,rowLength[iRow],rowLower[iRow],rowUpper[iRow]);
     2996             nPossibleZeroCost++;
     2997           } else if (value!=-COIN_DBL_MAX) {
     2998             if (model->logLevel() > 4)
     2999             printf("Singleton %d with objective in row with %d equal elements - rhs %g,%g\n",iColumn,rowLength[iRow],rowLower[iRow],rowUpper[iRow]);
     3000             nPossibleNonzeroCost++;
     3001           }
     3002         }
     3003       }
     3004       if (nPossibleZeroCost||nPossibleNonzeroCost)
     3005         printf("%d singletons with zero cost, %d with valid cost\n",
     3006                nPossibleZeroCost,nPossibleNonzeroCost);
     3007       // look for DW
     3008       int * blockStart = new int [2*(numberRows+numberColumns)+1+numberRows];
     3009       int * columnBlock = blockStart+numberRows;
     3010       int * nextColumn = columnBlock+numberColumns;
     3011       int * blockCount = nextColumn+numberColumns;
     3012       int * blockEls = blockCount+numberRows+1;
     3013       int direction[2]={-1,1};
     3014       int bestBreak=-1;
     3015       double bestValue=0.0;
     3016       int iPass=0;
     3017       int halfway=(numberRows+1)/2;
     3018       int firstMaster=-1;
     3019       int lastMaster=-2;
     3020       while (iPass<2) {
     3021         int increment=direction[iPass];
     3022         int start= increment>0 ? 0 : numberRows-1;
     3023         int stop=increment>0 ? numberRows : -1;
     3024         int numberBlocks=0;
     3025         int thisBestBreak=-1;
     3026         double thisBestValue=COIN_DBL_MAX;
     3027         int numberRowsDone=0;
     3028         int numberMarkedColumns=0;
     3029         int maximumBlockSize=0;
     3030         for (int i=0;i<numberRows+2*numberColumns;i++)
     3031           blockStart[i]=-1;
     3032         for (int i=0;i<numberRows+1;i++)
     3033           blockCount[i]=0;
     3034         for (int iRow=start;iRow!=stop;iRow+=increment) {
     3035           int iBlock = -1;
     3036           for (CoinBigIndex j=rowStart[iRow];j<rowStart[iRow]+rowLength[iRow];j++) {
     3037             int iColumn=column[j];
     3038             int whichColumnBlock=columnBlock[iColumn];
     3039             if (whichColumnBlock>=0) {
     3040               // column marked
     3041               if (iBlock<0) {
     3042                 // put row in that block
     3043                 iBlock=whichColumnBlock;
     3044               } else if (iBlock!=whichColumnBlock) {
     3045                 // merge
     3046                 blockCount[iBlock]+=blockCount[whichColumnBlock];
     3047                 blockCount[whichColumnBlock]=0;
     3048                 int jColumn=blockStart[whichColumnBlock];
     3049                 while (jColumn>=0) {
     3050                   columnBlock[jColumn]=iBlock;
     3051                   iColumn=jColumn;
     3052                   jColumn=nextColumn[jColumn];
     3053                 }
     3054                 nextColumn[iColumn]=blockStart[iBlock];
     3055                 blockStart[iBlock]=blockStart[whichColumnBlock];
     3056                 blockStart[whichColumnBlock]=-1;
     3057               }
     3058             }
     3059           }
     3060           int n=numberMarkedColumns;
     3061           if (iBlock<0) {
     3062             //new block
     3063             if (rowLength[iRow]) {
     3064               numberBlocks++;
     3065               iBlock=numberBlocks;
     3066               int jColumn=column[rowStart[iRow]];
     3067               columnBlock[jColumn]=iBlock;
     3068               blockStart[iBlock]=jColumn;
     3069               numberMarkedColumns++;
     3070               for (CoinBigIndex j=rowStart[iRow]+1;j<rowStart[iRow]+rowLength[iRow];j++) {
     3071                 int iColumn=column[j];
     3072                 columnBlock[iColumn]=iBlock;
     3073                 numberMarkedColumns++;
     3074                 nextColumn[jColumn]=iColumn;
     3075                 jColumn=iColumn;
     3076               }
     3077               blockCount[iBlock]=numberMarkedColumns-n;
     3078             } else {
     3079               // empty
     3080               iBlock=numberRows;
     3081             }
     3082           } else {
     3083             // put in existing block
     3084             int jColumn=blockStart[iBlock];
     3085             for (CoinBigIndex j=rowStart[iRow];j<rowStart[iRow]+rowLength[iRow];j++) {
     3086               int iColumn=column[j];
     3087               assert (columnBlock[iColumn]<0||columnBlock[iColumn]==iBlock);
     3088               if (columnBlock[iColumn]<0) {
     3089                 columnBlock[iColumn]=iBlock;
     3090                 numberMarkedColumns++;
     3091                 nextColumn[iColumn]=jColumn;
     3092                 jColumn=iColumn;
     3093               }
     3094             }
     3095             blockStart[iBlock]=jColumn;
     3096             blockCount[iBlock]+=numberMarkedColumns-n;
     3097           }
     3098           maximumBlockSize=CoinMax(maximumBlockSize,blockCount[iBlock]);
     3099           numberRowsDone++;
     3100           if (thisBestValue*numberRowsDone > maximumBlockSize&&numberRowsDone>halfway) {
     3101             thisBestBreak=iRow;
     3102             thisBestValue=static_cast<double>(maximumBlockSize)/
     3103               static_cast<double>(numberRowsDone);
     3104           }
     3105         }
     3106         if (thisBestBreak==stop)
     3107           thisBestValue=COIN_DBL_MAX;
     3108         iPass++;
     3109         if (iPass==1) {
     3110           bestBreak=thisBestBreak;
     3111           bestValue=thisBestValue;
     3112         } else {
     3113           if (bestValue<thisBestValue) {
     3114             firstMaster=0;
     3115             lastMaster=bestBreak;
     3116           } else {
     3117             firstMaster=thisBestBreak; // ? +1
     3118             lastMaster=numberRows;
     3119           }
     3120         }
     3121       }
     3122       if (firstMaster<lastMaster) {
     3123         printf("%d master rows %d <= < %d\n",lastMaster-firstMaster,
     3124                firstMaster,lastMaster);
     3125         for (int i=0;i<numberRows+2*numberColumns;i++)
     3126           blockStart[i]=-1;
     3127         for (int i=firstMaster;i<lastMaster;i++)
     3128           blockStart[i]=-2;
     3129         int firstRow=0;
     3130         int numberBlocks=-1;
     3131         while (true) {
     3132           for (;firstRow<numberRows;firstRow++) {
     3133             if (blockStart[firstRow]==-1)
     3134               break;
     3135           }
     3136           if (firstRow==numberRows)
     3137             break;
     3138           int nRows=0;
     3139           numberBlocks++;
     3140           int numberStack=1;
     3141           blockCount[0] = firstRow;
     3142           while (numberStack) {
     3143             int iRow=blockCount[--numberStack];
     3144             for (CoinBigIndex j=rowStart[iRow];j<rowStart[iRow]+rowLength[iRow];j++) {
     3145               int iColumn=column[j];
     3146               int iBlock=columnBlock[iColumn];
     3147               if (iBlock<0) {
     3148                 columnBlock[iColumn]=numberBlocks;
     3149                 for (CoinBigIndex k=columnStart[iColumn];
     3150                      k<columnStart[iColumn]+columnLength[iColumn];k++) {
     3151                   int jRow=row[k];
     3152                   int rowBlock=blockStart[jRow];
     3153                   if (rowBlock==-1) {
     3154                     nRows++;
     3155                     blockStart[jRow]=numberBlocks;
     3156                     blockCount[numberStack++]=jRow;
     3157                   }
     3158                 }
     3159               }
     3160             }
     3161           }
     3162           if (!nRows) {
     3163             // empty!!
     3164             numberBlocks--;
     3165           }
     3166           firstRow++;
     3167         }
     3168         // adjust
     3169         numberBlocks++;
     3170         for (int i=0;i<numberBlocks;i++) {
     3171           blockCount[i]=0;
     3172           nextColumn[i]=0;
     3173         }
     3174         int numberEmpty=0;
     3175         int numberMaster=0;
     3176         memset(blockEls,0,numberBlocks*sizeof(int));
     3177         for (int iRow = 0; iRow < numberRows; iRow++) {
     3178           int iBlock=blockStart[iRow];
     3179           if (iBlock>=0) {
     3180             blockCount[iBlock]++;
     3181             blockEls[iBlock]+=rowLength[iRow];
     3182           } else {
     3183             if (iBlock==-2)
     3184               numberMaster++;
     3185             else
     3186               numberEmpty++;
     3187           }
     3188         }
     3189         int numberEmptyColumns=0;
     3190         int numberMasterColumns=0;
     3191         for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
     3192           int iBlock=columnBlock[iColumn];
     3193           if (iBlock>=0) {
     3194             nextColumn[iBlock]++;
     3195           } else {
     3196             if (columnLength[iColumn])
     3197               numberMasterColumns++;
     3198             else
     3199               numberEmptyColumns++;
     3200           }
     3201         }
     3202         int largestRows=0;
     3203         int largestColumns=0;
     3204         for (int i=0;i<numberBlocks;i++) {
     3205           if (blockCount[i]+nextColumn[i]>largestRows+largestColumns) {
     3206             largestRows=blockCount[i];
     3207             largestColumns=nextColumn[i];
     3208           }
     3209         }
     3210         bool useful=true;
     3211         if (numberMaster>halfway||largestRows*3>numberRows)
     3212           useful=false;
     3213         printf("%s %d blocks (largest %d,%d), %d master rows (%d empty) out of %d, %d master columns (%d empty) out of %d\n",
     3214                useful ? "**Useful" : "NoGood",
     3215                numberBlocks,largestRows,largestColumns,numberMaster,numberEmpty,numberRows,
     3216                numberMasterColumns,numberEmptyColumns,numberColumns);
     3217         for (int i=0;i<numberBlocks;i++)
     3218           printf("Block %d has %d rows and %d columns (%d elements)\n",
     3219                  i,blockCount[i],nextColumn[i],blockEls[i]);
     3220         if (model->logLevel() == 17) {
     3221           int * whichRows=new int[numberRows+numberColumns];
     3222           int * whichColumns=whichRows+numberRows;
     3223           char name[20];
     3224           for (int iBlock=0;iBlock<numberBlocks;iBlock++) {
     3225             sprintf(name,"block%d.mps",iBlock);
     3226             int nRows=0;
     3227             for (int iRow=0;iRow<numberRows;iRow++) {
     3228               if (blockStart[iRow]==iBlock)
     3229                 whichRows[nRows++]=iRow;
     3230             }
     3231             int nColumns=0;
     3232             for (int iColumn=0;iColumn<numberColumns;iColumn++) {
     3233               if (columnBlock[iColumn]==iBlock)
     3234                 whichColumns[nColumns++]=iColumn;
     3235             }
     3236             ClpSimplex subset(model,nRows,whichRows,nColumns,whichColumns);
     3237             subset.writeMps(name,0,1);
     3238           }
     3239           delete [] whichRows;
     3240         }
     3241       }
     3242       delete [] blockStart;
     3243     }
    27093244     for (iRow = 0; iRow < numberRows; iRow++) {
    27103245          int length = rowLength[iRow];
     
    34013936  1.11.00 November 5 2009 (Guy Fawkes) - OSL factorization and better ordering
    34023937 */
     3938#ifdef CILK_TEST
     3939// -*- C++ -*-
     3940
     3941/*
     3942 * cilk-for.cilk
     3943 *
     3944 * Copyright (c) 2007-2008 Cilk Arts, Inc.  55 Cambridge Street,
     3945 * Burlington, MA 01803.  Patents pending.  All rights reserved. You may
     3946 * freely use the sample code to guide development of your own works,
     3947 * provided that you reproduce this notice in any works you make that
     3948 * use the sample code.  This sample code is provided "AS IS" without
     3949 * warranty of any kind, either express or implied, including but not
     3950 * limited to any implied warranty of non-infringement, merchantability
     3951 * or fitness for a particular purpose.  In no event shall Cilk Arts,
     3952 * Inc. be liable for any direct, indirect, special, or consequential
     3953 * damages, or any other damages whatsoever, for any use of or reliance
     3954 * on this sample code, including, without limitation, any lost
     3955 * opportunity, lost profits, business interruption, loss of programs or
     3956 * data, even if expressly advised of or otherwise aware of the
     3957 * possibility of such damages, whether in an action of contract,
     3958 * negligence, tort, or otherwise.
     3959 *
     3960 * This file demonstrates a Cilk++ for loop
     3961 */
     3962
     3963#include <cilk/cilk.h>
     3964//#include <cilk/cilkview.h>
     3965#include <cilk/reducer_max.h>
     3966#include <cstdlib>
     3967#include <iostream>
     3968
     3969// cilk_for granularity.
     3970#define CILK_FOR_GRAINSIZE 128
     3971
     3972double dowork(double i)
     3973{
     3974    // Waste time:
     3975    int j;
     3976    double k = i;
     3977    for (j = 0; j < 50000; ++j) {
     3978        k += k / ((j + 1) * (k + 1));
     3979    }
     3980
     3981    return k;
     3982}
     3983static void doSomeWork(double * a,int low, int high)
     3984{
     3985  if (high-low>300) {
     3986    int mid=(high+low)>>1;
     3987    cilk_spawn doSomeWork(a,low,mid);
     3988    doSomeWork(a,mid,high);
     3989    cilk_sync;
     3990  } else {
     3991    for(int i = low; i < high; ++i) {
     3992      a[i] = dowork(a[i]);
     3993    }
     3994  }
     3995}
     3996
     3997void cilkTest()
     3998{
     3999    unsigned int n = 10000;
     4000    //cilk::cilkview cv;
     4001
     4002
     4003    double* a = new double[n];
     4004
     4005    for(unsigned int i = 0; i < n; i++) {
     4006        // Populate A
     4007        a[i] = (double) ((i * i) % 1024 + 512) / 512;
     4008    }
     4009
     4010    std::cout << "Iterating over " << n << " integers" << std::endl;
     4011
     4012    //cv.start();
     4013#if 1
     4014    //#pragma cilk_grainsize=CILK_FOR_GRAINSIZE
     4015    cilk_for(unsigned int i = 0; i < n; ++i) {
     4016        a[i] = dowork(a[i]);
     4017    }
     4018#else
     4019    doSomeWork(a,0,n);
     4020#endif
     4021    int * which =new int[n];
     4022    unsigned int n2=n>>1;
     4023    for (int i=0;i<n2;i++)
     4024      which[i]=n-2*i;
     4025    cilk::reducer_max_index<int,double> maximumIndex(-1,0.0);
     4026    cilk_for(unsigned int i = 0; i < n2; ++i) {
     4027      int iWhich=which[i];
     4028      maximumIndex.calc_max(iWhich,a[iWhich]);
     4029    }
     4030    int bestIndex=maximumIndex.get_index();
     4031    int bestIndex2=-1;
     4032    double largest=0.0;
     4033    cilk_for(unsigned int i = 0; i < n2; ++i) {
     4034      int iWhich=which[i];
     4035      if (a[iWhich]>largest) {
     4036        bestIndex2=iWhich;
     4037        largest=a[iWhich];
     4038      }
     4039    }
     4040    assert (bestIndex==bestIndex2);
     4041    //cv.stop();
     4042    //cv.dump("cilk-for-results", false);
     4043
     4044    //std::cout << cv.accumulated_milliseconds() / 1000.f << " seconds" << std::endl;
     4045
     4046    exit(0);
     4047}
     4048#endif
  • trunk/Clp/src/ClpMessage.cpp

    r1817 r1878  
    2222     {CLP_SIMPLEX_STATUS, 6, 1, "%d  Obj %g%? Primal inf %g (%d)%? Dual inf %g (%d)%? w.o. free dual inf (%d)"},
    2323     {CLP_DUAL_BOUNDS, 25, 3, "Looking optimal checking bounds with %g"},
     24#if ABC_NORMAL_DEBUG>1
     25     {CLP_SIMPLEX_ACCURACY, 60, 1, "Primal error %g, dual error %g"},
     26#else
    2427     {CLP_SIMPLEX_ACCURACY, 60, 3, "Primal error %g, dual error %g"},
     28#endif
    2529     {CLP_SIMPLEX_BADFACTOR, 7, 2, "Singular factorization of basis - status %d"},
    2630     {CLP_SIMPLEX_BOUNDTIGHTEN, 8, 3, "Bounds were tightened %d times"},
  • trunk/Clp/src/ClpModel.hpp

    r1825 r1878  
    356356     }
    357357#ifndef CLP_NO_STD
    358      inline std::string problemName() const {
     358     inline const std::string & problemName() const {
    359359          return strParam_[ClpProbName];
    360360     }
     
    11761176       shift by 65336 is 3 all same, 1 all except col bounds
    11771177     */
     1178#define ROW_COLUMN_COUNTS_SAME 1
     1179#define MATRIX_SAME 2
     1180#define MATRIX_JUST_ROWS_ADDED 4
     1181#define MATRIX_JUST_COLUMNS_ADDED 8
     1182#define ROW_LOWER_SAME 16
     1183#define ROW_UPPER_SAME 32
     1184#define OBJECTIVE_SAME 64
     1185#define COLUMN_LOWER_SAME 128
     1186#define COLUMN_UPPER_SAME 256
     1187#define BASIS_SAME 512
     1188#define ALL_SAME 65339
     1189#define ALL_SAME_EXCEPT_COLUMN_BOUNDS 65337
    11781190     unsigned int whatsChanged_;
    11791191     /// Status of problem
  • trunk/Clp/src/ClpPackedMatrix.cpp

    r1865 r1878  
    40874087                         printf("Out of range %d %lld %d %g\n", iColumn, j, row[j], elementByColumn[j]);
    40884088#endif
    4089                          delete [] mark;
    40904089                         return false;
    40914090                    }
     
    48364835          if (matrix_->isColOrdered() && numberOther > matrix_->getNumCols())
    48374836               matrix_->setDimensions(-1, numberOther);
    4838           if (!matrix_->isColOrdered() || numberOther >= 0 || matrix_->getExtraGap() || matrix_->hasGaps()) {
     4837          if (!matrix_->isColOrdered() || numberOther >= 0 || matrix_->getExtraGap()) {
    48394838               numberErrors = matrix_->appendRows(number, starts, index, element, numberOther);
    48404839          } else {
  • trunk/Clp/src/ClpPresolve.cpp

    r1834 r1878  
    1313
    1414#include "CoinHelperFunctions.hpp"
     15#if CLP_HAS_ABC
     16#include "CoinAbcCommon.hpp"
     17#endif
    1518
    1619#include "CoinPackedMatrix.hpp"
     
    5558     nrows_(0),
    5659     nelems_(0),
     60#ifdef ABC_INHERIT
     61     numberPasses_(20),
     62#else
    5763     numberPasses_(5),
     64#endif
    5865     substitution_(3),
    5966#ifndef CLP_NO_STD
     
    274281#ifndef CLP_NO_STD
    275282     }
    276 #endif
     283#endif 
    277284     // put back duals
    278285     CoinMemcpyN(prob.rowduals_,        nrows_, originalModel_->dualRowSolution());
     
    444451}
    445452#endif
     453static int tightenDoubletons2(CoinPresolveMatrix * prob)
     454{
     455  // column-major representation
     456  const int ncols = prob->ncols_ ;
     457  const CoinBigIndex *const mcstrt = prob->mcstrt_ ;
     458  const int *const hincol = prob->hincol_ ;
     459  const int *const hrow = prob->hrow_ ;
     460  double * colels = prob->colels_ ;
     461  double * cost = prob->cost_ ;
     462
     463  // column type, bounds, solution, and status
     464  const unsigned char *const integerType = prob->integerType_ ;
     465  double * clo = prob->clo_ ;
     466  double * cup = prob->cup_ ;
     467  // row-major representation
     468  //const int nrows = prob->nrows_ ;
     469  const CoinBigIndex *const mrstrt = prob->mrstrt_ ;
     470  const int *const hinrow = prob->hinrow_ ;
     471  const int *const hcol = prob->hcol_ ;
     472  double * rowels = prob->rowels_ ;
     473
     474  // row bounds
     475  double *const rlo = prob->rlo_ ;
     476  double *const rup = prob->rup_ ;
     477
     478  // tolerances
     479  //const double ekkinf2 = PRESOLVE_SMALL_INF ;
     480  //const double ekkinf = ekkinf2*1.0e8 ;
     481  //const double ztolcbarj = prob->ztoldj_ ;
     482  //const CoinRelFltEq relEq(prob->ztolzb_) ;
     483  int numberChanged=0;
     484  double bound[2];
     485  double alpha[2]={0.0,0.0};
     486  double offset=0.0;
     487
     488  for (int icol=0;icol<ncols;icol++) {
     489    if (hincol[icol]==2) {
     490      CoinBigIndex start=mcstrt[icol];
     491      int row0 = hrow[start];
     492      if (hinrow[row0]!=2)
     493        continue;
     494      int row1 = hrow[start+1];
     495      if (hinrow[row1]!=2)
     496        continue;
     497      double element0 = colels[start];
     498      double rowUpper0=rup[row0];
     499      bool swapSigns0=false;
     500      if (rlo[row0]>-1.0e30) {
     501        if (rup[row0]>1.0e30) {
     502          swapSigns0=true;
     503          rowUpper0=-rlo[row0];
     504          element0=-element0;
     505        } else {
     506          // range or equality
     507          continue;
     508        }
     509      } else if (rup[row0]>1.0e30) {
     510        // free
     511        continue;
     512      }
     513#if 0
     514      // skip here for speed
     515      // skip if no cost (should be able to get rid of)
     516      if (!cost[icol]) {
     517        printf("should be able to get rid of %d with no cost\n",icol);
     518        continue;
     519      }
     520      // skip if negative cost for now
     521      if (cost[icol]<0.0) {
     522        printf("code for negative cost\n");
     523        continue;
     524      }
     525#endif
     526      double element1 = colels[start+1];
     527      double rowUpper1=rup[row1];
     528      bool swapSigns1=false;
     529      if (rlo[row1]>-1.0e30) {
     530        if (rup[row1]>1.0e30) {
     531          swapSigns1=true;
     532          rowUpper1=-rlo[row1];
     533          element1=-element1;
     534        } else {
     535          // range or equality
     536          continue;
     537        }
     538      } else if (rup[row1]>1.0e30) {
     539        // free
     540        continue;
     541      }
     542      double lowerX=clo[icol];
     543      double upperX=cup[icol];
     544      int otherCol=-1;
     545      CoinBigIndex startRow=mrstrt[row0];
     546      for (CoinBigIndex j=startRow;j<startRow+2;j++) {
     547        int jcol=hcol[j];
     548        if (jcol!=icol) {
     549          alpha[0]=swapSigns0 ? -rowels[j] :rowels[j];
     550          otherCol=jcol;
     551        }
     552      }
     553      startRow=mrstrt[row1];
     554      bool possible=true;
     555      for (CoinBigIndex j=startRow;j<startRow+2;j++) {
     556        int jcol=hcol[j];
     557        if (jcol!=icol) {
     558          if (jcol==otherCol) {
     559            alpha[1]=swapSigns1 ? -rowels[j] :rowels[j];
     560          } else {
     561            possible=false;
     562          }
     563        }
     564      }
     565      if (possible) {
     566        // skip if no cost (should be able to get rid of)
     567        if (!cost[icol]) {
     568          printf("should be able to get rid of %d with no cost\n",icol);
     569          continue;
     570        }
     571        // skip if negative cost for now
     572        if (cost[icol]<0.0) {
     573          printf("code for negative cost\n");
     574          continue;
     575        }
     576        bound[0]=clo[otherCol];
     577        bound[1]=cup[otherCol];
     578        double lowestLowest=COIN_DBL_MAX;
     579        double highestLowest=-COIN_DBL_MAX;
     580        double lowestHighest=COIN_DBL_MAX;
     581        double highestHighest=-COIN_DBL_MAX;
     582        int binding0=0;
     583        int binding1=0;
     584        for (int k=0;k<2;k++) {
     585          bool infLow0=false;
     586          bool infLow1=false;
     587          double sum0=0.0;
     588          double sum1=0.0;
     589          double value=bound[k];
     590          if (fabs(value)<1.0e30) {
     591            sum0+=alpha[0]*value;
     592            sum1+=alpha[1]*value;
     593          } else {
     594            if (alpha[0]>0.0) {
     595              if (value<0.0)
     596                infLow0 =true;
     597            } else if (alpha[0]<0.0) {
     598              if (value>0.0)
     599                infLow0 =true;
     600            }
     601            if (alpha[1]>0.0) {
     602              if (value<0.0)
     603                infLow1 =true;
     604            } else if (alpha[1]<0.0) {
     605              if (value>0.0)
     606                infLow1 =true;
     607            }
     608          }
     609          /* Got sums
     610           */
     611          double thisLowest0=-COIN_DBL_MAX;
     612          double thisHighest0=COIN_DBL_MAX;
     613          if (element0>0.0) {
     614            // upper bound unless inf&2 !=0
     615            if (!infLow0)
     616              thisHighest0 = (rowUpper0-sum0)/element0;
     617          } else {
     618            // lower bound unless inf&2 !=0
     619            if (!infLow0)
     620              thisLowest0 = (rowUpper0-sum0)/element0;
     621          }
     622          double thisLowest1=-COIN_DBL_MAX;
     623          double thisHighest1=COIN_DBL_MAX;
     624          if (element1>0.0) {
     625            // upper bound unless inf&2 !=0
     626            if (!infLow1)
     627              thisHighest1 = (rowUpper1-sum1)/element1;
     628          } else {
     629            // lower bound unless inf&2 !=0
     630            if (!infLow1)
     631              thisLowest1 = (rowUpper1-sum1)/element1;
     632          }
     633          if (thisLowest0>thisLowest1+1.0e-12) {
     634            if (thisLowest0>lowerX+1.0e-12)
     635              binding0|= 1<<k;
     636          } else if (thisLowest1>thisLowest0+1.0e-12) {
     637            if (thisLowest1>lowerX+1.0e-12)
     638              binding1|= 1<<k;
     639            thisLowest0=thisLowest1;
     640          }
     641          if (thisHighest0<thisHighest1-1.0e-12) {
     642            if (thisHighest0<upperX-1.0e-12)
     643              binding0|= 1<<k;
     644          } else if (thisHighest1<thisHighest0-1.0e-12) {
     645            if (thisHighest1<upperX-1.0e-12)
     646              binding1|= 1<<k;
     647            thisHighest0=thisHighest1;
     648          }
     649          lowestLowest=CoinMin(lowestLowest,thisLowest0);
     650          highestHighest=CoinMax(highestHighest,thisHighest0);
     651          lowestHighest=CoinMin(lowestHighest,thisHighest0);
     652          highestLowest=CoinMax(highestLowest,thisLowest0);
     653        }
     654        // see if any good
     655        //#define PRINT_VALUES
     656        if (!binding0||!binding1) {
     657          printf("Row redundant for column %d\n",icol);
     658        } else {
     659#ifdef PRINT_VALUES
     660          printf("Column %d bounds %g,%g lowest %g,%g highest %g,%g\n",
     661                 icol,lowerX,upperX,lowestLowest,lowestHighest,
     662                 highestLowest,highestHighest);
     663#endif
     664          // if integer adjust
     665          if (integerType[icol]) {
     666            lowestLowest=ceil(lowestLowest-1.0e-5);
     667            highestLowest=ceil(highestLowest-1.0e-5);
     668            lowestHighest=floor(lowestHighest+1.0e-5);
     669            highestHighest=floor(highestHighest+1.0e-5);
     670          }
     671          // if costed may be able to adjust
     672          if (cost[icol]>=0.0) {
     673            if (highestLowest<upperX&&highestLowest>=lowerX&&highestHighest<1.0e30) {
     674              highestHighest=CoinMin(highestHighest,highestLowest);
     675            }
     676          }
     677          if (cost[icol]<=0.0) {
     678            if (lowestHighest>lowerX&&lowestHighest<=upperX&&lowestHighest>-1.0e30) {
     679              lowestLowest=CoinMax(lowestLowest,lowestHighest);
     680            }
     681          }
     682#if 1
     683          if (lowestLowest>lowerX+1.0e-8) {
     684#ifdef PRINT_VALUES
     685            printf("Can increase lower bound on %d from %g to %g\n",
     686                   icol,lowerX,lowestLowest);
     687#endif
     688            lowerX=lowestLowest;
     689          }
     690          if (highestHighest<upperX-1.0e-8) {
     691#ifdef PRINT_VALUES
     692            printf("Can decrease upper bound on %d from %g to %g\n",
     693                   icol,upperX,highestHighest);
     694#endif
     695            upperX=highestHighest;
     696           
     697          }
     698#endif
     699          // see if we can move costs
     700          double xValue;
     701          double yValue0;
     702          double yValue1;
     703          double newLower=COIN_DBL_MAX;
     704          double newUpper=-COIN_DBL_MAX;
     705          double ranges0[2];
     706          double ranges1[2];
     707          double costEqual;
     708          double slope[2];
     709          assert (binding0+binding1==3);
     710          // get where equal
     711          xValue=(rowUpper0*element1-rowUpper1*element0)/(alpha[0]*element1-alpha[1]*element0);
     712          yValue0=(rowUpper0-xValue*alpha[0])/element0;
     713          yValue1=(rowUpper1-xValue*alpha[1])/element1;
     714          newLower=CoinMin(newLower,CoinMax(yValue0,yValue1));
     715          newUpper=CoinMax(newUpper,CoinMax(yValue0,yValue1));
     716          double xValueEqual=xValue;
     717          double yValueEqual=yValue0;
     718          costEqual = xValue*cost[otherCol]+yValueEqual*cost[icol];
     719          if (binding0==1) {
     720            ranges0[0]=bound[0];
     721            ranges0[1]=yValue0;
     722            ranges1[0]=yValue0;
     723            ranges1[1]=bound[1];
     724            // take x 1.0 down
     725            double x=xValue-1.0;
     726            double y=(rowUpper0-x*alpha[0])/element0;
     727            double costTotal = x*cost[otherCol]+y*cost[icol];
     728            slope[0] = costEqual-costTotal;
     729            // take x 1.0 up
     730            x=xValue+1.0;
     731            y=(rowUpper1-x*alpha[1])/element0;
     732            costTotal = x*cost[otherCol]+y*cost[icol];
     733            slope[1] = costTotal-costEqual;
     734          } else {
     735            ranges1[0]=bound[0];
     736            ranges1[1]=yValue0;
     737            ranges0[0]=yValue0;
     738            ranges0[1]=bound[1];
     739            // take x 1.0 down
     740            double x=xValue-1.0;
     741            double y=(rowUpper1-x*alpha[1])/element0;
     742            double costTotal = x*cost[otherCol]+y*cost[icol];
     743            slope[1] = costEqual-costTotal;
     744            // take x 1.0 up
     745            x=xValue+1.0;
     746            y=(rowUpper0-x*alpha[0])/element0;
     747            costTotal = x*cost[otherCol]+y*cost[icol];
     748            slope[0] = costTotal-costEqual;
     749          }
     750#ifdef PRINT_VALUES
     751          printf("equal value of %d is %g, value of %d is max(%g,%g) - %g\n",
     752                 otherCol,xValue,icol,yValue0,yValue1,CoinMax(yValue0,yValue1));
     753          printf("Cost at equality %g for constraint 0 ranges %g -> %g slope %g for constraint 1 ranges %g -> %g slope %g\n",
     754                 costEqual,ranges0[0],ranges0[1],slope[0],ranges1[0],ranges1[1],slope[1]);
     755#endif
     756          xValue=bound[0];
     757          yValue0=(rowUpper0-xValue*alpha[0])/element0;
     758          yValue1=(rowUpper1-xValue*alpha[1])/element1;
     759#ifdef PRINT_VALUES
     760          printf("value of %d is %g, value of %d is max(%g,%g) - %g\n",
     761                 otherCol,xValue,icol,yValue0,yValue1,CoinMax(yValue0,yValue1));
     762#endif
     763          newLower=CoinMin(newLower,CoinMax(yValue0,yValue1));
     764          // cost>0 so will be at lower
     765          //double yValueAtBound0=newLower;
     766          newUpper=CoinMax(newUpper,CoinMax(yValue0,yValue1));
     767          xValue=bound[1];
     768          yValue0=(rowUpper0-xValue*alpha[0])/element0;
     769          yValue1=(rowUpper1-xValue*alpha[1])/element1;
     770#ifdef PRINT_VALUES
     771          printf("value of %d is %g, value of %d is max(%g,%g) - %g\n",
     772                 otherCol,xValue,icol,yValue0,yValue1,CoinMax(yValue0,yValue1));
     773#endif
     774          newLower=CoinMin(newLower,CoinMax(yValue0,yValue1));
     775          // cost>0 so will be at lower
     776          //double yValueAtBound1=newLower;
     777          newUpper=CoinMax(newUpper,CoinMax(yValue0,yValue1));
     778          lowerX=CoinMax(lowerX,newLower-1.0e-12*fabs(newLower));
     779          upperX=CoinMin(upperX,newUpper+1.0e-12*fabs(newUpper));
     780          // Now make duplicate row
     781          // keep row 0 so need to adjust costs so same
     782#ifdef PRINT_VALUES
     783          printf("Costs for x %g,%g,%g are %g,%g,%g\n",
     784                 xValueEqual-1.0,xValueEqual,xValueEqual+1.0,
     785                 costEqual-slope[0],costEqual,costEqual+slope[1]);
     786#endif
     787          double costOther=cost[otherCol]+slope[1];
     788          double costThis=cost[icol]+slope[1]*(element0/alpha[0]);
     789          xValue=xValueEqual;
     790          yValue0=CoinMax((rowUpper0-xValue*alpha[0])/element0,lowerX);
     791          double thisOffset=costEqual-(costOther*xValue+costThis*yValue0);
     792          offset += thisOffset;
     793#ifdef PRINT_VALUES
     794          printf("new cost at equal %g\n",costOther*xValue+costThis*yValue0+thisOffset);
     795#endif
     796          xValue=xValueEqual-1.0;
     797          yValue0=CoinMax((rowUpper0-xValue*alpha[0])/element0,lowerX);
     798#ifdef PRINT_VALUES
     799          printf("new cost at -1 %g\n",costOther*xValue+costThis*yValue0+thisOffset);
     800#endif
     801          assert(fabs((costOther*xValue+costThis*yValue0+thisOffset)-(costEqual-slope[0]))<1.0e-5);
     802          xValue=xValueEqual+1.0;
     803          yValue0=CoinMax((rowUpper0-xValue*alpha[0])/element0,lowerX);
     804#ifdef PRINT_VALUES
     805          printf("new cost at +1 %g\n",costOther*xValue+costThis*yValue0+thisOffset);
     806#endif
     807          assert(fabs((costOther*xValue+costThis*yValue0+thisOffset)-(costEqual+slope[1]))<1.0e-5);
     808          numberChanged++;
     809          //      continue;
     810          cost[otherCol] = costOther;
     811          cost[icol] = costThis;
     812          clo[icol]=lowerX;
     813          cup[icol]=upperX;
     814          int startCol[2];
     815          int endCol[2];
     816          startCol[0]=mcstrt[icol];
     817          endCol[0]=startCol[0]+2;
     818          startCol[1]=mcstrt[otherCol];
     819          endCol[1]=startCol[1]+hincol[otherCol];
     820          double values[2]={0.0,0.0};
     821          for (int k=0;k<2;k++) {
     822            for (CoinBigIndex i=startCol[k];i<endCol[k];i++) {
     823              if (hrow[i]==row0)
     824                values[k]=colels[i];
     825            }
     826            for (CoinBigIndex i=startCol[k];i<endCol[k];i++) {
     827              if (hrow[i]==row1)
     828                colels[i]=values[k];
     829            }
     830          }
     831          for (CoinBigIndex i=mrstrt[row1];i<mrstrt[row1]+2;i++) {
     832            if (hcol[i]==icol)
     833              rowels[i]=values[0];
     834            else
     835              rowels[i]=values[1];
     836          }
     837        }
     838      }
     839    }
     840  }
     841#if ABC_NORMAL_DEBUG>0
     842  if (offset)
     843    printf("Cost offset %g\n",offset);
     844#endif
     845  return numberChanged;
     846}
    446847//#define COIN_PRESOLVE_BUG
    447848#ifdef COIN_PRESOLVE_BUG
     
    557958          }
    558959
     960            // transfer costs (may want to do it in OsiPresolve)
     961            // need a transfer back at end of postsolve transferCosts(prob);
    559962
    560963          int iLoop;
     
    569972               paction_ = dupcol_action::presolve(prob, paction_);
    570973          }
     974#ifdef ABC_INHERIT
     975          if (doTwoxTwo()) {
     976            possibleSkip;
     977            paction_ = twoxtwo_action::presolve(prob, paction_);
     978          }
     979#endif
    571980          if (duprow) {
    572981            possibleSkip;
    573                paction_ = duprow_action::presolve(prob, paction_);
     982            int nTightened=tightenDoubletons2(prob);
     983            if (nTightened)
     984              printf("%d doubletons tightened\n",nTightened);
     985            paction_ = duprow_action::presolve(prob, paction_);
    574986          }
    575987          if (doGubrow()) {
     
    583995          int lastDropped = 0;
    584996          prob->pass_ = 0;
     997#ifdef ABC_INHERIT
     998          int numberRowsStart=nrows_-prob->countEmptyRows();
     999          int numberColumnsStart=ncols_-prob->countEmptyCols();
     1000          int numberRowsLeft=numberRowsStart;
     1001          int numberColumnsLeft=numberColumnsStart;
     1002          bool lastPassWasGood=true;
     1003#if ABC_NORMAL_DEBUG
     1004          printf("Original rows,columns %d,%d starting first pass with %d,%d\n",
     1005                 nrows_,ncols_,numberRowsLeft,numberColumnsLeft);
     1006#endif
     1007#endif
    5851008          if (numberPasses_<=5)
    5861009              prob->presolveOptions_ |= 0x10000; // say more lightweight
     
    8671290               if (paction_ == paction0 || stopLoop)
    8681291                    break;
    869 
     1292#ifdef ABC_INHERIT
     1293               // see whether to stop anyway
     1294               int numberRowsNow=nrows_-prob->countEmptyRows();
     1295               int numberColumnsNow=ncols_-prob->countEmptyCols();
     1296#if ABC_NORMAL_DEBUG
     1297               printf("Original rows,columns %d,%d - last %d,%d end of pass %d has %d,%d\n",
     1298                      nrows_,ncols_,numberRowsLeft,numberColumnsLeft,iLoop+1,numberRowsNow,
     1299                      numberColumnsNow);
     1300#endif
     1301               int rowsDeleted=numberRowsLeft-numberRowsNow;
     1302               int columnsDeleted=numberColumnsLeft-numberColumnsNow;
     1303               if (iLoop>15) {
     1304                 if (rowsDeleted*100<numberRowsStart&&
     1305                     columnsDeleted*100<numberColumnsStart)
     1306                   break;
     1307                 lastPassWasGood=true;
     1308               } else if (rowsDeleted*100<numberRowsStart&&rowsDeleted<500&&
     1309                          columnsDeleted*100<numberColumnsStart&&columnsDeleted<500) {
     1310                 if (!lastPassWasGood)
     1311                   break;
     1312                 else
     1313                   lastPassWasGood=false;
     1314               } else {
     1315                 lastPassWasGood=true;
     1316               }
     1317               numberRowsLeft=numberRowsNow;
     1318               numberColumnsLeft=numberColumnsNow;
     1319#endif
    8701320          }
    8711321     }
     
    9431393               }
    9441394          }
     1395     }
     1396     if (prob.maxmin_<0) {
     1397       //for (int i=0;i<presolvedModel_->numberRows();i++)
     1398       //prob.rowduals_[i]=-prob.rowduals_[i];
     1399       for (int i=0;i<ncols_;i++) {
     1400         prob.cost_[i]=-prob.cost_[i];
     1401         //prob.rcosts_[i]=-prob.rcosts_[i];
     1402       }
     1403       prob.maxmin_=1.0;
    9451404     }
    9461405     const CoinPresolveAction *paction = paction_;
     
    12451704     mcstrt_[0] = 0;
    12461705     ClpDisjointCopyN(m->getVectorLengths(), ncols_,  hincol_);
     1706     if (si->getObjSense() < 0.0) {
     1707       for (int i=0;i<ncols_;i++)
     1708         cost_[i]=-cost_[i];
     1709       maxmin_=1.0;
     1710     }
    12471711     for (icol = 0; icol < ncols_; icol++) {
    12481712          CoinBigIndex j;
     
    14161880#endif
    14171881{
     1882     if (si->getObjSense() < 0.0) {
     1883       for (int i=0;i<ncols_;i++)
     1884         cost_[i]=-cost_[i];
     1885       dobias_=-dobias_;
     1886     }
    14181887     si->loadProblem(ncols_, nrows_, mcstrt_, hrow_, colels_, hincol_,
    14191888                     clo_, cup_, cost_, rlo_, rup_);
     
    14361905#endif
    14371906     si->setDblParam(ClpObjOffset, originalOffset_ - dobias_);
     1907     if (si->getObjSense() < 0.0) {
     1908       // put back
     1909       for (int i=0;i<ncols_;i++)
     1910         cost_[i]=-cost_[i];
     1911       dobias_=-dobias_;
     1912       maxmin_=-1.0;
     1913     }
    14381914
    14391915}
     
    19042380               if (rowObjective_) {
    19052381                    int iRow;
     2382#ifndef NDEBUG
    19062383                    int k = -1;
     2384#endif
    19072385                    int nObj = 0;
    19082386                    for (iRow = 0; iRow < nrowsNow; iRow++) {
    19092387                         int kRow = originalRow_[iRow];
     2388#ifndef NDEBUG
    19102389                         assert (kRow > k);
    19112390                         k = kRow;
     2391#endif
    19122392                         rowObjective_[iRow] = rowObjective_[kRow];
    19132393                         if (rowObjective_[iRow])
  • trunk/Clp/src/ClpPresolve.hpp

    r1780 r1878  
    168168          else presolveActions_ |= 1024;
    169169     }
     170     /// Whether we want to do twoxtwo part of presolve
     171     inline bool doTwoxTwo() const {
     172          return (presolveActions_ & 2048) == 0;
     173     }
     174     inline void setDoTwoxtwo(bool doTwoxTwo) {
     175          if (doTwoxTwo) presolveActions_  &= ~2048;
     176          else presolveActions_ |= 2048;
     177     }
    170178     /// Set whole group
    171179     inline int presolveActions() const {
  • trunk/Clp/src/ClpPrimalColumnSteepest.cpp

    r1732 r1878  
    591591          }
    592592     }
     593#endif
     594#if 0
     595     for (int i=0;i<numberRows;i++)
     596       printf("row %d weight %g infeas %g\n",i,weights_[i+numberColumns],infeas[i+numberColumns]);
     597     for (int i=0;i<numberColumns;i++)
     598       printf("column %d weight %g infeas %g\n",i,weights_[i],infeas[i]);
    593599#endif
    594600     return bestSequence;
  • trunk/Clp/src/ClpSimplex.cpp

    r1872 r1878  
    3232#include "CoinLpIO.hpp"
    3333#include <cfloat>
     34#if CLP_HAS_ABC
     35#include "CoinAbcCommon.hpp"
     36#endif
    3437
    3538#include <string>
     
    122125     perturbationArray_(NULL),
    123126     baseModel_(NULL)
     127#ifdef ABC_INHERIT
     128     ,abcSimplex_(NULL),
     129     abcState_(0)
     130#endif
    124131{
    125132     int i;
     
    233240     perturbationArray_(NULL),
    234241     baseModel_(NULL)
     242#ifdef ABC_INHERIT
     243     ,abcSimplex_(NULL),
     244     abcState_(0)
     245#endif
    235246{
    236247     int i;
     
    381392     perturbationArray_(NULL),
    382393     baseModel_(NULL)
     394#ifdef ABC_INHERIT
     395     ,abcSimplex_(NULL),
     396     abcState_(rhs->abcState_)
     397#endif
    383398{
    384399     int i;
     
    17421757     }
    17431758#    endif
     1759#if 0 //ndef NDEBUG
     1760  // Make sure everything is clean
     1761  double sumOutside=0.0;
     1762  int numberOutside=0;
     1763  //double sumOutsideLarge=0.0;
     1764  int numberOutsideLarge=0;
     1765  double sumInside=0.0;
     1766  int numberInside=0;
     1767  //double sumInsideLarge=0.0;
     1768  int numberInsideLarge=0;
     1769  int numberTotal=numberRows_+numberColumns_;
     1770  for (int iSequence = 0; iSequence < numberTotal; iSequence++) {
     1771    if(getStatus(iSequence) == isFixed) {
     1772      // double check fixed
     1773      assert (upper_[iSequence] == lower_[iSequence]);
     1774      assert (fabs(solution_[iSequence]-lower_[iSequence])<primalTolerance_);
     1775    } else if (getStatus(iSequence) == isFree) {
     1776      assert (upper_[iSequence] == COIN_DBL_MAX && lower_[iSequence]==-COIN_DBL_MAX);
     1777    } else if (getStatus(iSequence) == atLowerBound) {
     1778      assert (fabs(solution_[iSequence]-lower_[iSequence])<1000.0*primalTolerance_);
     1779      if (solution_[iSequence]<lower_[iSequence]) {
     1780        numberOutside++;
     1781        sumOutside-=solution_[iSequence]-lower_[iSequence];
     1782        if (solution_[iSequence]<lower_[iSequence]-primalTolerance_)
     1783          numberOutsideLarge++;
     1784      } else if (solution_[iSequence]>lower_[iSequence]) {
     1785        numberInside++;
     1786        sumInside+=solution_[iSequence]-lower_[iSequence];
     1787        if (solution_[iSequence]>lower_[iSequence]+primalTolerance_)
     1788          numberInsideLarge++;
     1789      }
     1790    } else if (getStatus(iSequence) == atUpperBound) {
     1791      assert (fabs(solution_[iSequence]-upper_[iSequence])<1000.0*primalTolerance_);
     1792      if (solution_[iSequence]>upper_[iSequence]) {
     1793        numberOutside++;
     1794        sumOutside+=solution_[iSequence]-upper_[iSequence];
     1795        if (solution_[iSequence]>upper_[iSequence]+primalTolerance_)
     1796          numberOutsideLarge++;
     1797      } else if (solution_[iSequence]<upper_[iSequence]) {
     1798        numberInside++;
     1799        sumInside-=solution_[iSequence]-upper_[iSequence];
     1800        if (solution_[iSequence]<upper_[iSequence]-primalTolerance_)
     1801          numberInsideLarge++;
     1802      }
     1803    } else if (getStatus(iSequence) == superBasic) {
     1804      //assert (!valuesPass);
     1805    }
     1806  }
     1807  if (numberInside+numberOutside)
     1808    printf("%d outside summing to %g (%d large), %d inside summing to %g (%d large)\n",
     1809           numberOutside,sumOutside,numberOutsideLarge,
     1810           numberInside,sumInside,numberInsideLarge);
     1811#endif
    17441812     int status = factorization_->factorize(this, solveType, valuesPass);
    17451813     if (status) {
     
    20232091          return 1;
    20242092     } else if (numberIterations_ > 1000 + 10 * (numberRows_ + (numberColumns_ >> 2)) && matrix_->type()<15) {
    2025           // A bit worried problem may be cycling - lets factorize at random intervals for a short period
    2026           int numberTooManyIterations = numberIterations_ - 10 * (numberRows_ + (numberColumns_ >> 2));
    20272093          double random = randomNumberGenerator_.randomDouble();
    2028           int window = numberTooManyIterations%5000;
    2029           if (window<2*maximumPivots)
    2030             random = 0.2*random+0.8; // randomly re-factorize but not too soon
    2031           else
    2032             random=1.0; // switch off if not in window of opportunity
    20332094          int maxNumber = (forceFactorization_ < 0) ? maximumPivots : CoinMin(forceFactorization_, maximumPivots);
    20342095          if (factorization_->pivots() >= random * maxNumber) {
     
    21302191     perturbationArray_(NULL),
    21312192     baseModel_(NULL)
     2193#ifdef ABC_INHERIT
     2194     ,abcSimplex_(NULL),
     2195     abcState_(0)
     2196#endif
    21322197{
    21332198     int i;
     
    22342299     perturbationArray_(NULL),
    22352300     baseModel_(NULL)
     2301#ifdef ABC_INHERIT
     2302     ,abcSimplex_(NULL),
     2303     abcState_(0)
     2304#endif
    22362305{
    22372306     int i;
     
    24102479     allowedInfeasibility_ = rhs.allowedInfeasibility_;
    24112480     automaticScale_ = rhs.automaticScale_;
     2481#ifdef ABC_INHERIT
     2482     abcSimplex_ = NULL;
     2483     abcState_ = rhs.abcState_;
     2484#endif
    24122485     maximumPerturbationSize_ = rhs.maximumPerturbationSize_;
    24132486     if (maximumPerturbationSize_ && maximumPerturbationSize_ >= 2 * numberColumns_) {
     
    31043177                    maximumInternalColumns_ = numberColumns_;
    31053178               }
     3179               //maximumRows_=CoinMax(maximumInternalRows_,maximumRows_);
     3180               //maximumColumns_=CoinMax(maximumInternalColumns_,maximumColumns_);
    31063181               assert(maximumInternalRows_ == maximumRows_);
    31073182               assert(maximumInternalColumns_ == maximumColumns_);
     
    32033278          inverseRowScale_ = NULL;
    32043279          inverseColumnScale_ = NULL;
     3280          if (scalingFlag_ > 0 &&(specialOptions_ & 65536) != 0&&
     3281              rowScale_&&rowScale_==savedRowScale_)
     3282            rowScale_=NULL;
    32053283          if (scalingFlag_ > 0 && !rowScale_) {
    32063284               if ((specialOptions_ & 65536) != 0) {
     
    38063884          for (int i = 0; i < numberRows2; i++) {
    38073885               int iSequence = pivotVariable_[i];
    3808                if (iSequence < numberRows_ + numberColumns_ &&
    3809                          getStatus(iSequence) != basic) {
     3886               if (iSequence<0||(iSequence < numberRows_ + numberColumns_ &&
     3887                                 getStatus(iSequence) != basic)) {
    38103888                    keepPivots = false;
    38113889                    break;
     
    38903968                         columnArray_[iColumn] = new CoinIndexedVector();
    38913969                    }
    3892                     if (!iColumn)
    3893                          columnArray_[iColumn]->reserve(numberColumns_);
    3894                     else
    3895                          columnArray_[iColumn]->reserve(CoinMax(numberRows2, numberColumns_));
     3970                    columnArray_[iColumn]->reserve(numberColumns_+numberRows2);
    38963971               }
    38973972          } else {
     
    74447519     }
    74457520     if (!numberBasic || pivot == 3) {
     7521#if 0
     7522          int nFree=0;
     7523          for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
     7524            if (columnLower_[iColumn] < -1.0e20 && columnUpper_[iColumn] > 1.0e20)
     7525              nFree++;
     7526          }
     7527          if (nFree>numberColumns_/10)
     7528            pivot=3;
     7529#endif
    74467530          if (pivot == 3) {
    7447                // Just throw all free variables in basis
    7448                for (int iRow = 0; iRow < numberRows_; iRow++) {
    7449                     if (fabs(rowLower_[iRow]) < fabs(rowUpper_[iRow]))
    7450                          setRowStatus(iRow, atLowerBound);
    7451                     else
    7452                          setRowStatus(iRow, atUpperBound);
    7453                }
    7454                for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
    7455                     if (columnLower_[iColumn] < -1.0e20 && columnUpper_[iColumn] > 1.0e20)
    7456                          setColumnStatus(iColumn, basic);
    7457                }
    7458                return 0;
     7531            // Get column copy
     7532            CoinPackedMatrix * columnCopy = matrix();
     7533            const int * row = columnCopy->getIndices();
     7534            const CoinBigIndex * columnStart = columnCopy->getVectorStarts();
     7535            const int * columnLength = columnCopy->getVectorLengths();
     7536            const double * element = columnCopy->getElements();
     7537            int nFree=0;
     7538            for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
     7539              if (columnLower_[iColumn] < -1.0e20 && columnUpper_[iColumn] > 1.0e20) {
     7540                // find triangular row
     7541                int kRow=-1;
     7542                double largest=0.0;
     7543                double largestOther=0.0;
     7544                for (CoinBigIndex j=columnStart[iColumn];
     7545                     j<columnStart[iColumn]+columnLength[iColumn];j++) {
     7546                  int iSequence=row[j]+numberColumns_;
     7547                  if (!flagged(iSequence)) {
     7548                    if (fabs(element[j])>largest) {
     7549                      kRow=row[j];
     7550                      largest=fabs(element[j]);
     7551                    }
     7552                  } else {
     7553                    if (fabs(element[j])>largestOther) {
     7554                      largestOther=fabs(element[j]);
     7555                    }
     7556                  }
     7557                }
     7558                if (kRow>=0&&largest*2.5>=largestOther) {
     7559                  nFree++;
     7560                  setColumnStatus(iColumn, basic);
     7561                  if (fabs(rowLower_[kRow]) < fabs(rowUpper_[kRow]))
     7562                    setRowStatus(kRow, atLowerBound);
     7563                  else
     7564                    setRowStatus(kRow, atUpperBound);
     7565                  for (CoinBigIndex j=columnStart[iColumn];
     7566                       j<columnStart[iColumn]+columnLength[iColumn];j++) {
     7567                    int iSequence=row[j]+numberColumns_;
     7568                    setFlagged(iSequence);
     7569                  }
     7570                }
     7571              }
     7572            }
     7573            if (nFree) {
     7574              for (int i=0;i<numberRows_;i++)
     7575                clearFlagged(i);
     7576              printf("%d free variables put in basis\n",nFree);
     7577              return 0;
     7578            }
    74597579          }
    74607580          // all slack
     
    78787998                              << numberBad
    78797999                              << CoinMessageEol;
     8000                    // could do elsewhere and could even if something done
     8001                    // see if some MUST be in basis
     8002                    if (!numberIn) {
     8003                    }
    78808004                    returnCode =  -1;
    78818005               }
     
    82708394                    if (algorithm_ > 0 && (objective_->type() < 2 || !objective_->activated())) {
    82718395#ifndef FEB_TRY
    8272                          static_cast<ClpSimplexPrimal *> (this)->perturb(0);
     8396                      //static_cast<ClpSimplexPrimal *> (this)->perturb(0);
    82738397#endif
    82748398                    } else if (algorithm_ < 0) {
     
    83618485          } else {
    83628486               // using previous factorization - we assume fine
    8363                if ((moreSpecialOptions_ & 8) == 0 || !(whatsChanged_ & 64)) {
    8364                     // but we need to say not optimal (!(whatsChanged_ & 64) means that objective has changed)
     8487               if ((moreSpecialOptions_ & 8) == 0) {
     8488                    // but we need to say not optimal
    83658489                    numberPrimalInfeasibilities_ = 1;
    83668490                    numberDualInfeasibilities_ = 1;
     
    97329856                    // dual should be zero
    97339857                    if (fabs(dualValue) > 10.0 * dualTolerance) {
    9734                          sumDualInfeasibilities_ -= dualValue + dualTolerance_;
     9858                      sumDualInfeasibilities_ += fabs(dualValue) - dualTolerance_;
    97359859                         numberDualInfeasibilities_ ++;
    97369860                    }
     
    993410058          // User did not touch preset
    993510059          const int cutoff1 = 10000;
    9936           const int cutoff2 = 100000;
    993710060          const int base = 75;
    993810061          const int freq0 = 50;
     10062          int frequency;
     10063#ifndef ABC_INHERIT
     10064          const int cutoff2 = 100000;
    993910065          const int freq1 = 200;
    994010066          const int freq2 = 400;
    994110067          const int maximum = 1000;
    9942           int frequency;
    994310068          if (numberRows_ < cutoff1)
    994410069               frequency = base + numberRows_ / freq0;
     
    994710072          else
    994810073               frequency = base + cutoff1 / freq0 + (cutoff2 - cutoff1) / freq1 + (numberRows_ - cutoff2) / freq2;
     10074#else
     10075          const int freq1 = 150;
     10076          const int maximum = 10000;
     10077          if (numberRows_ < cutoff1)
     10078               frequency = base + numberRows_ / freq0;
     10079          else
     10080               frequency = base + cutoff1 / freq0 + (numberRows_ - cutoff1) / freq1;
     10081#endif
    994910082          setFactorizationFrequency(CoinMin(maximum, frequency));
    995010083     }
  • trunk/Clp/src/ClpSimplex.hpp

    r1868 r1878  
    2828class ClpDisasterHandler;
    2929class ClpConstraint;
    30 
     30#if CLP_HAS_ABC
     31#include "AbcCommon.hpp"
     32class AbcTolerancesEtc;
     33class AbcSimplex;
     34#include "CoinAbcCommon.hpp"
     35#endif
    3136/** This solves LPs using the simplex method
    3237
     
    116121         Only to be used with mini constructor */
    117122     void originalModel(ClpSimplex * miniModel);
     123  inline int abcState() const
     124  { return abcState_;}
     125  inline void setAbcState(int state)
     126  { abcState_=state;}
     127#ifdef ABC_INHERIT
     128  inline AbcSimplex * abcSimplex() const
     129  { return abcSimplex_;}
     130  inline void setAbcSimplex(AbcSimplex * simplex)
     131  { abcSimplex_=simplex;}
     132  /// Returns 0 if dual can be skipped
     133  int doAbcDual();
     134  /// Returns 0 if primal can be skipped
     135  int doAbcPrimal(int ifValuesPass);
     136#endif
    118137     /** Array persistence flag
    119138         If 0 then as now (delete/new)
     
    213232     int loadNonLinear(void * info, int & numberConstraints,
    214233                       ClpConstraint ** & constraints);
     234#ifdef ABC_INHERIT
     235  /// Loads tolerances etc
     236  void loadTolerancesEtc(const AbcTolerancesEtc & data);
     237  /// Unloads tolerances etc
     238  void unloadTolerancesEtc(AbcTolerancesEtc & data);
     239#endif
    215240     //@}
    216241
     
    276301     /// Solve using structure of model and maybe in parallel
    277302     int solve(CoinStructuredModel * model);
     303#ifdef ABC_INHERIT
     304  /** solvetype 0 for dual, 1 for primal
     305      startup 1 for values pass
     306      interrupt whether to pass across interrupt handler
     307  */
     308  void dealWithAbc(int solveType,int startUp,bool interrupt=false);
     309#endif
    278310     /** This loads a model from a CoinStructuredModel object - returns number of errors.
    279311         If originalOrder then keep to order stored in blocks,
     
    11471179         4096 bit - try more for complete fathoming
    11481180         8192 bit - don't even think of using primal if user asks for dual (and vv)
     1181         16384 bit - in initialSolve so be more flexible
     1182         debug
     1183         32768 bit - do dual in netlibd
     1184         65536 (*3) initial stateDualColumn
    11491185     */
    11501186     inline int moreSpecialOptions() const {
     
    11661202         4096 bit - try more for complete fathoming
    11671203         8192 bit - don't even think of using primal if user asks for dual (and vv)
     1204         16384 bit - in initialSolve so be more flexible
    11681205     */
    11691206     inline void setMoreSpecialOptions(int value) {
     
    15991636     /// For dealing with all issues of cycling etc
    16001637     ClpSimplexProgress progress_;
     1638#ifdef ABC_INHERIT
     1639  AbcSimplex * abcSimplex_;
     1640#define CLP_ABC_WANTED 1
     1641#define CLP_ABC_WANTED_PARALLEL 2
     1642#define CLP_ABC_FULL_DONE 8
     1643#define CLP_ABC_BEEN_FEASIBLE 65536
     1644  // bits 256,512,1024 for crash
     1645#endif
     1646  int abcState_;
    16011647public:
    16021648     /// Spare int array for passing information [0]!=0 switches on
  • trunk/Clp/src/ClpSimplexDual.cpp

    r1868 r1878  
    569569ClpSimplexDual::dual(int ifValuesPass, int startFinishOptions)
    570570{
     571  //handler_->setLogLevel(63);
     572  //yprintf("STARTing dual %d rows\n",numberRows_);
    571573     bestObjectiveValue_ = -COIN_DBL_MAX;
    572574     algorithm_ = -1;
     
    26282630     }
    26292631#endif
    2630      double freeAlpha = 0.0;
     2632     //double freeAlpha = 0.0;
    26312633     if (alreadyChosen < 0) {
    26322634          // first see if any free variables and put them in basis
     
    26782680               if (chosenRow >= 0) {
    26792681                    pivotRow_ = chosenRow;
    2680                     freeAlpha = work[chosenRow];
     2682                    //freeAlpha = work[chosenRow];
    26812683               }
    26822684               rowArray_[1]->clear();
     
    51955197               } else if (largestPrimalError_ > 1.0e5) {
    51965198                    {
    5197                          int iBigB = -1;
     5199                      //int iBigB = -1;
    51985200                         double bigB = 0.0;
    5199                          int iBigN = -1;
     5201                         //int iBigN = -1;
    52005202                         double bigN = 0.0;
    52015203                         for (int i = 0; i < numberRows_ + numberColumns_; i++) {
     
    52045206                                   if (value > bigB) {
    52055207                                        bigB = value;
    5206                                         iBigB = i;
     5208                                        //iBigB = i;
    52075209                                   }
    52085210                              } else {
    52095211                                   if (value > bigN) {
    52105212                                        bigN = value;
    5211                                         iBigN = i;
     5213                                        //iBigN = i;
    52125214                                   }
    52135215                              }
     
    54285430     return modified;
    54295431}
     5432#if ABC_NORMAL_DEBUG>0
     5433//#define PERT_STATISTICS
     5434#endif
     5435#ifdef PERT_STATISTICS
     5436static void breakdown(const char * name, int numberLook, const double * region)
     5437{
     5438     double range[] = {
     5439          -COIN_DBL_MAX,
     5440          -1.0e15, -1.0e11, -1.0e8, -1.0e5, -1.0e4, -1.0e3, -1.0e2, -1.0e1,
     5441          -1.0,
     5442          -1.0e-1, -1.0e-2, -1.0e-3, -1.0e-4, -1.0e-5, -1.0e-8, -1.0e-11, -1.0e-15,
     5443          0.0,
     5444          1.0e-15, 1.0e-11, 1.0e-8, 1.0e-5, 1.0e-4, 1.0e-3, 1.0e-2, 1.0e-1,
     5445          1.0,
     5446          1.0e1, 1.0e2, 1.0e3, 1.0e4, 1.0e5, 1.0e8, 1.0e11, 1.0e15,
     5447          COIN_DBL_MAX
     5448     };
     5449     int nRanges = static_cast<int> (sizeof(range) / sizeof(double));
     5450     int * number = new int[nRanges];
     5451     memset(number, 0, nRanges * sizeof(int));
     5452     int * numberExact = new int[nRanges];
     5453     memset(numberExact, 0, nRanges * sizeof(int));
     5454     int i;
     5455     for ( i = 0; i < numberLook; i++) {
     5456          double value = region[i];
     5457          for (int j = 0; j < nRanges; j++) {
     5458               if (value == range[j]) {
     5459                    numberExact[j]++;
     5460                    break;
     5461               } else if (value < range[j]) {
     5462                    number[j]++;
     5463                    break;
     5464               }
     5465          }
     5466     }
     5467     printf("\n%s has %d entries\n", name, numberLook);
     5468     for (i = 0; i < nRanges; i++) {
     5469          if (number[i])
     5470               printf("%d between %g and %g", number[i], range[i-1], range[i]);
     5471          if (numberExact[i]) {
     5472               if (number[i])
     5473                    printf(", ");
     5474               printf("%d exactly at %g", numberExact[i], range[i]);
     5475          }
     5476          if (number[i] + numberExact[i])
     5477               printf("\n");
     5478     }
     5479     delete [] number;
     5480     delete [] numberExact;
     5481}
     5482#endif
    54305483// Perturbs problem
    54315484int
     
    56805733     double weight[] = {1.0e-4, 1.0e-2, 5.0e-1, 1.0, 2.0, 5.0, 10.0, 20.0, 30.0, 40.0, 100.0};
    56815734     //double weight[]={1.0e-4,1.0e-2,5.0e-1,1.0,20.0,50.0,100.0,120.0,130.0,140.0,200.0};
    5682      double extraWeight = 10.0;
     5735     //double extraWeight = 10.0;
    56835736     // Scale back if wanted
    56845737     double weight2[] = {1.0e-4, 1.0e-2, 5.0e-1, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0};
    56855738     if (constantPerturbation < 99.0 * dualTolerance_) {
    56865739          perturbation *= 0.1;
    5687           extraWeight = 0.5;
     5740          //extraWeight = 0.5;
    56885741          memcpy(weight, weight2, sizeof(weight2));
    56895742     }
     
    58085861     //CoinZeroN(cost_+nTotal,nTotal);
    58095862     // say perturbed
     5863#ifdef PERT_STATISTICS
     5864  {
     5865    double averageCost = 0.0;
     5866    int numberNonZero = 0;
     5867    double * COIN_RESTRICT sort = new double[numberColumns_];
     5868    for (int i = 0; i < numberColumns_; i++) {
     5869      double value = fabs(cost_[i]);
     5870      sort[i] = value;
     5871      averageCost += value;
     5872      if (value)
     5873        numberNonZero++;
     5874    }
     5875    if (numberNonZero)
     5876      averageCost /= static_cast<double> (numberNonZero);
     5877    else
     5878      averageCost = 1.0;
     5879    std::sort(sort, sort + numberColumns_);
     5880    int number = 1;
     5881    double last = sort[0];
     5882    for (int i = 1; i < numberColumns_; i++) {
     5883      if (last != sort[i])
     5884        number++;
     5885    last = sort[i];
     5886    }
     5887    printf("nnz %d percent %d", number, (number * 100) / numberColumns_);
     5888    delete [] sort;
     5889    breakdown("Objective", numberColumns_+numberRows_, cost_);
     5890  }
     5891#endif
    58105892     perturbation_ = 101;
    58115893     return 0;
  • trunk/Clp/src/ClpSimplexOther.cpp

    r1877 r1878  
    2929#include <stdio.h>
    3030#include <iostream>
    31 #ifdef HAS_CILK
    32 #include <cilk/cilk.h>
    33 #else
    34 #define cilk_for for
    35 #define cilk_spawn
    36 #define cilk_sync
    37 #endif
    3831#ifdef INT_IS_8
    3932#define COIN_ANY_BITS_PER_INT 64
     
    121114               matrix_->transposeTimes(this, -1.0,
    122115                                       rowArray_[0], columnArray_[1], columnArray_[0]);
    123 #if COIN_FAC_NEW
     116#ifdef COIN_FAC_NEW
    124117               assert (!rowArray_[0]->packedMode());
    125118#endif
     
    138131               } else {
    139132                    int number = rowArray_[0]->getNumElements();
    140 #if COIN_FAC_NEW
     133#ifdef COIN_FAC_NEW
    141134                    const int * index = rowArray_[0]->getIndices();
    142135#endif
     
    807800     bool changed = false;
    808801     int numberChanged = 0;
     802     int numberFreeColumnsInPrimal=0;
    809803     int iColumn;
    810804     // check if we need to change bounds to rows
    811805     for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
    812           if (columnUpper_[iColumn] < 1.0e20 &&
    813                     columnLower_[iColumn] > -1.0e20) {
     806       if (columnUpper_[iColumn] < 1.0e20) {
     807         if (columnLower_[iColumn] > -1.0e20) {
    814808               changed = true;
    815809               numberChanged++;
    816810          }
     811       } else if (columnLower_[iColumn] < -1.0e20) {
     812         numberFreeColumnsInPrimal++;
     813       }
    817814     }
    818815     int iRow;
    819816     int numberExtraRows = 0;
     817     int numberFreeColumnsInDual=0;
    820818     if (numberChanged <= fractionColumnRanges * numberColumns_) {
    821819          for (iRow = 0; iRow < numberRows_; iRow++) {
     
    824822                    if (rowUpper_[iRow] != rowLower_[iRow])
    825823                         numberExtraRows++;
     824                    else
     825                      numberFreeColumnsInDual++;
    826826               }
    827827          }
     
    831831          return NULL;
    832832     }
     833     printf("would have %d free columns in primal, %d in dual\n",
     834            numberFreeColumnsInPrimal,numberFreeColumnsInDual);
     835     if (4*(numberFreeColumnsInDual-numberFreeColumnsInPrimal)>
     836         numberColumns_&&fractionRowRanges<1.0)
     837       return NULL; //dangerous (well anyway in dual)
    833838     if (changed) {
    834839          ClpSimplex * model3 = new ClpSimplex(*model2);
     
    978983     modelDual->setSpecialOptions(model2->specialOptions());
    979984     modelDual->setMoreSpecialOptions(model2->moreSpecialOptions());
     985     modelDual->setMaximumIterations(model2->maximumIterations());
     986     modelDual->setFactorizationFrequency(model2->factorizationFrequency());
     987     modelDual->setLogLevel(model2->logLevel());
    980988     delete [] fromRowsLower;
    981989     delete [] fromRowsUpper;
     
    11661174                    setRowStatus(iRow, atUpperBound);
    11671175               } else {
    1168                     assert (dualDj[iRow] < 1.0e-5);
     1176                    // might be stopped assert (dualDj[iRow] < 1.0e-5);
    11691177                    rowActivity_[iRow] = rowUpper_[iRow] + dualDj[iRow];
    11701178               }
     
    11751183               } else {
    11761184                    rowActivity_[iRow] = rowLower_[iRow] + dualDj[iRow];
    1177                     assert (dualDj[iRow] > -1.0e-5);
     1185                    // might be stopped assert (dualDj[iRow] > -1.0e-5);
    11781186               }
    11791187          } else {
     
    11911199                    //     dualSol[kExtraRow],dualDj[iRow],dualDj[kExtraRow]);
    11921200                    if (status == basic) {
    1193                          assert (statusL != basic);
     1201                         // might be stopped assert (statusL != basic);
    11941202                         rowActivity_[iRow] = rowUpper_[iRow];
    11951203                         setRowStatus(iRow, atUpperBound);
     
    12011209                    } else {
    12021210                         rowActivity_[iRow] = rowLower_[iRow] - dualDj[iRow];
    1203                          assert (dualDj[iRow] < 1.0e-5);
     1211                         // might be stopped assert (dualDj[iRow] < 1.0e-5);
    12041212                         // row basic
    12051213                         //setRowStatus(iRow,basic);
     
    33923400   
    33933401    // exit if victory declared
    3394     if ((problemStatus_ >= 0 && startingTheta>=endingTheta-1.0e-7) ||
    3395         (problemStatus_ == 0 && startingTheta>=CoinMin(endingTheta-1.0e-7,1.0e49) )) {
    3396       startingTheta=endingTheta;
     3402    if (problemStatus_ >= 0 && startingTheta>=endingTheta-1.0e-7 )
    33973403      break;
    3398     }
    33993404   
    34003405    // test for maximum iterations
     
    76887693#endif
    76897694  double inverseCleanup = (cleanUp>0.0) ? 1.0/cleanUp : 0.0;
    7690   /*
    7691     Rules for cleaning up rhs
    7692     if cleanup != 0.0 just use cleanup
    7693     otherwise if problemStatus_==0 (i.e. was optimal) use rowActivity
    7694     then go to integer value if close
    7695    */
    76967695  //#define PRINT_DUP
    76977696#ifdef PRINT_DUP
     
    78417840            break;
    78427841          } else if (fabs(rup2-rlo2)<=tolerance) {
    7843             // equal
    7844             if (!inverseCleanup) {
    7845               if (!problemStatus_) {
    7846                 // choose one closer to rowActivity
    7847                 double rval2=rowActivity_[iThis];
    7848                 if (fabs(rup2-rval2)<fabs(rlo2-rval2))
    7849                   rlo2=rup2;
    7850               } else {
    7851                 // no solution - closer to zero
    7852                 if (fabs(rup2)<fabs(rlo2))
    7853                   rlo2=rup2;
    7854               }
    7855               // see if close to integer value
    7856               double value2 = floor(rlo2+0.5);
    7857               if (fabs(rlo2-value2)<1.0e-9)
    7858                 rlo2=value2;
     7842            // equal - choose closer to zero
     7843            if (fabs(rup2)<fabs(rlo2))
     7844              rlo2=rup2;
     7845            else
    78597846              rup2=rlo2;
    7860             } else {
    7861               // we will be cleaning up later - choose closer to zero
    7862               if (fabs(rup2)<fabs(rlo2))
    7863                 rlo2=rup2;
    7864               else
    7865                 rup2=rlo2;
    7866             }
    78677847#if 0
    78687848          if (rowLength[iThis]<4)
     
    87288708        double ratio = fabs(value1/value2);
    87298709        if (ratio<0.001||ratio>1000.0) {
    8730           if (fabs(value1)>fabs(value2)) {
     8710          if (fabs(value1)<fabs(value2)) {
    87318711            swap=true;
    87328712          }
     
    98799859        int length=columnLengthX[iColumn];
    98809860        CoinBigIndex start=columnStartX[iColumn];
    9881         int newLength=length+((jRowLower<0||jRowUpper<0) ? 1 : 2);
    98829861        int nextColumn=forward[iColumn];
    98839862        CoinBigIndex startNext=columnStartX[nextColumn];
    9884         if (start+newLength>startNext) {
     9863        if (start+length==startNext) {
    98859864          // need more
    98869865          moveAround(numberColumns_,numberElementsOriginal,
    9887                      iColumn,newLength,
     9866                     iColumn,length+(jRowLower<0||jRowUpper<0) ? 1 : 2,
    98889867                     forward,backward,columnStartX,columnLengthX,
    98899868                     rowX,elementX);
  • trunk/Clp/src/ClpSimplexPrimal.cpp

    r1868 r1878  
    497497                         }
    498498                         //#define FEB_TRY
    499 #ifdef FEB_TRY
     499#if 1 //def FEB_TRY
    500500                         if (perturbation_ < 100)
    501501                              perturb(0);
     
    740740#ifdef CLP_USER_DRIVEN
    741741               // If large number of pivots trap later?
    742                if (problemStatus_==-1 && numberPivots<1000) {
     742               if (problemStatus_==-1 && numberPivots<1000000) {
    743743                 int status = eventHandler_->event(ClpEventHandler::noCandidateInPrimal);
    744744                 if (status>=0&&status<10) {
     
    12801280                         goToDual = unPerturb(); // stop any further perturbation
    12811281                         if (nonLinearCost_->sumInfeasibilities() > 1.0e-1)
    1282                               goToDual = false;
     1282                           goToDual = false;
    12831283                         nonLinearCost_->checkInfeasibilities(primalTolerance_);
    12841284                         numberDualInfeasibilities_ = 1; // carry on
     
    15241524     if (numberIterations_ > lastBadIteration_ + 100)
    15251525          moreSpecialOptions_ &= ~16; // clear check accuracy flag
     1526     if ((moreSpecialOptions_ & 256) != 0)
     1527       goToDual=false;
    15261528     if (goToDual || (numberIterations_ > 1000 && largestPrimalError_ > 1.0e6
    15271529                      && largestDualError_ > 1.0e6)) {
     
    17401742               }
    17411743               double value;
    1742                assert (oldValue >= -tolerance);
     1744               //assert (oldValue >= -10.0*tolerance);
    17431745               if (possible) {
    17441746                    value = oldValue - upperTheta * alpha;
     
    23012303          cleanStatus(); // make sure status okay
    23022304     // Make sure feasible bounds
    2303      if (nonLinearCost_)
    2304           nonLinearCost_->feasibleBounds();
     2305     if (nonLinearCost_) {
     2306       nonLinearCost_->checkInfeasibilities();
     2307       //nonLinearCost_->feasibleBounds();
     2308     }
    23052309     // look at element range
    23062310     double smallestNegative;
     
    24652469          maximumFraction *= 0.1;
    24662470     }
     2471     if (savePerturbation==51) {
     2472       perturbation = CoinMin(0.1,perturbation);
     2473       maximumFraction *=0.1;
     2474     }
    24672475     if (number != numberRows_)
    24682476          type = 1;
     
    24962504                         double value = maximumFraction * (difference + bias);
    24972505                         value = CoinMin(value, 0.1);
     2506                         value = CoinMax(value,primalTolerance_);
    24982507#ifndef SAVE_PERT
    24992508                         value *= randomNumberGenerator_.randomDouble();
  • trunk/Clp/src/ClpSolve.cpp

    r1794 r1878  
    3737#include "ClpMessage.hpp"
    3838#include "CoinTime.hpp"
     39#if CLP_HAS_ABC
     40#include "CoinAbcCommon.hpp"
     41#endif
     42#ifdef ABC_INHERIT
     43#include "AbcSimplex.hpp"
     44#include "AbcSimplexFactorization.hpp"
     45#endif
     46double zz_slack_value=0.0;
    3947
    4048#include "ClpPresolve.hpp"
     
    5563#include "CoinPackedMatrix.hpp"
    5664#include "CoinMpsIO.hpp"
    57 
    5865//#############################################################################
    5966
     
    340347#include "CoinSignal.hpp"
    341348static ClpSimplex * currentModel = NULL;
     349#ifdef ABC_INHERIT
     350static AbcSimplex * currentAbcModel = NULL;
     351#endif
    342352
    343353extern "C" {
     
    350360          if (currentModel != NULL)
    351361               currentModel->setMaximumIterations(0); // stop at next iterations
     362#ifdef ABC_INHERIT
     363          if (currentAbcModel != NULL)
     364               currentAbcModel->setMaximumIterations(0); // stop at next iterations
     365#endif
    352366#ifndef SLIM_CLP
    353367          if (currentModel2 != NULL)
     
    357371     }
    358372}
    359 
     373#if ABC_INSTRUMENT>1
     374int abcPricing[20];
     375int abcPricingDense[20];
     376static int trueNumberRows;
     377static int numberTypes;
     378#define MAX_TYPES 25
     379#define MAX_COUNT 20
     380#define MAX_FRACTION 101
     381static char * types[MAX_TYPES];
     382static double counts[MAX_TYPES][MAX_COUNT];
     383static double countsFraction[MAX_TYPES][MAX_FRACTION];
     384static double * currentCounts;
     385static double * currentCountsFraction;
     386static int currentType;
     387static double workMultiplier[MAX_TYPES];
     388static double work[MAX_TYPES];
     389static double currentWork;
     390static double otherWork[MAX_TYPES];
     391static int timesCalled[MAX_TYPES];
     392static int timesStarted[MAX_TYPES];
     393static int fractionDivider;
     394void instrument_initialize(int numberRows)
     395{
     396  trueNumberRows=numberRows;
     397  numberTypes=0;
     398  memset(counts,0,sizeof(counts));
     399  currentCounts=NULL;
     400  memset(countsFraction,0,sizeof(countsFraction));
     401  currentCountsFraction=NULL;
     402  memset(workMultiplier,0,sizeof(workMultiplier));
     403  memset(work,0,sizeof(work));
     404  memset(otherWork,0,sizeof(otherWork));
     405  memset(timesCalled,0,sizeof(timesCalled));
     406  memset(timesStarted,0,sizeof(timesStarted));
     407  currentType=-1;
     408  fractionDivider=(numberRows+MAX_FRACTION-2)/(MAX_FRACTION-1);
     409}
     410void instrument_start(const char * type,int numberRowsEtc)
     411{
     412  if (currentType>=0)
     413    instrument_end();
     414  currentType=-1;
     415  currentWork=0.0;
     416  for (int i=0;i<numberTypes;i++) {
     417    if (!strcmp(types[i],type)) {
     418      currentType=i;
     419      break;
     420    }
     421  }
     422  if (currentType==-1) {
     423    assert (numberTypes<MAX_TYPES);
     424    currentType=numberTypes;
     425    types[numberTypes++]=strdup(type);
     426  }
     427  currentCounts = &counts[currentType][0];
     428  currentCountsFraction = &countsFraction[currentType][0];
     429  timesStarted[currentType]++;
     430  assert (trueNumberRows);
     431  workMultiplier[currentType]+=static_cast<double>(numberRowsEtc)/static_cast<double>(trueNumberRows);
     432}
     433void instrument_add(int count)
     434{
     435  assert (currentType>=0);
     436  currentWork+=count;
     437  timesCalled[currentType]++;
     438  if (count<MAX_COUNT-1)
     439    currentCounts[count]++;
     440  else
     441    currentCounts[MAX_COUNT-1]++;
     442  assert(count/fractionDivider>=0&&count/fractionDivider<MAX_FRACTION);
     443  currentCountsFraction[count/fractionDivider]++;
     444}
     445void instrument_do(const char * type,double count)
     446{
     447  int iType=-1;
     448  for (int i=0;i<numberTypes;i++) {
     449    if (!strcmp(types[i],type)) {
     450      iType=i;
     451      break;
     452    }
     453  }
     454  if (iType==-1) {
     455    assert (numberTypes<MAX_TYPES);
     456    iType=numberTypes;
     457    types[numberTypes++]=strdup(type);
     458  }
     459  timesStarted[iType]++;
     460  otherWork[iType]+=count;
     461}
     462void instrument_end()
     463{
     464  work[currentType]+=currentWork;
     465  currentType=-1;
     466}
     467void instrument_end_and_adjust(double factor)
     468{
     469  work[currentType]+=currentWork*factor;
     470  currentType=-1;
     471}
     472void instrument_print()
     473{
     474  for (int iType=0;iType<numberTypes;iType++) {
     475    currentCounts = &counts[iType][0];
     476    currentCountsFraction = &countsFraction[iType][0];
     477    if (!otherWork[iType]) {
     478      printf("%s started %d times, used %d times, work %g (average length %.1f) multiplier %g\n",
     479             types[iType],timesStarted[iType],timesCalled[iType],
     480             work[iType],work[iType]/(timesCalled[iType]+1.0e-100),workMultiplier[iType]/(timesStarted[iType]+1.0e-100));
     481      int n=0;
     482      for (int i=0;i<MAX_COUNT-1;i++) {
     483        if (currentCounts[i]) {
     484          if (n==5) {
     485            n=0;
     486            printf("\n");
     487          }
     488          n++;
     489          printf("(%d els,%.0f times) ",i,currentCounts[i]);
     490        }
     491      }
     492      if (currentCounts[MAX_COUNT-1]) {
     493        if (n==5) {
     494          n=0;
     495          printf("\n");
     496        }
     497        n++;
     498        printf("(>=%d els,%.0f times) ",MAX_COUNT-1,currentCounts[MAX_COUNT-1]);
     499      }
     500      printf("\n");
     501      int largestFraction;
     502      int nBig=0;
     503      for (largestFraction=MAX_FRACTION-1;largestFraction>=10;largestFraction--) {
     504        double count = currentCountsFraction[largestFraction];
     505        if (count&&largestFraction>10)
     506          nBig++;
     507        if (nBig>4)
     508          break;
     509      }
     510      int chunk=(largestFraction+5)/10;
     511      int lo=0;
     512      for (int iChunk=0;iChunk<largestFraction;iChunk+=chunk) {
     513        int hi=CoinMin(lo+chunk*fractionDivider,trueNumberRows);
     514        double sum=0.0;
     515        for (int i=iChunk;i<CoinMin(iChunk+chunk,MAX_FRACTION);i++)
     516          sum += currentCountsFraction[i];
     517        if (sum)
     518          printf("(%d-%d %.0f) ",lo,hi,sum);
     519        lo=hi;
     520      }
     521      for (int i=lo/fractionDivider;i<MAX_FRACTION;i++) {
     522        if (currentCountsFraction[i])
     523          printf("(%d %.0f) ",i*fractionDivider,currentCountsFraction[i]);
     524      }
     525      printf("\n");
     526    } else {
     527      printf("%s started %d times, used %d times, work %g multiplier %g other work %g\n",
     528             types[iType],timesStarted[iType],timesCalled[iType],
     529             work[iType],workMultiplier[iType],otherWork[iType]);
     530    }
     531    free(types[iType]);
     532  }
     533}
     534#endif
     535#if ABC_PARALLEL==2
     536#ifndef FAKE_CILK
     537int number_cilk_workers=0;
     538#include <cilk/cilk_api.h>
     539#endif
     540#endif
     541#ifdef ABC_INHERIT
     542void
     543ClpSimplex::dealWithAbc(int solveType, int startUp,
     544                        bool interrupt)
     545{
     546  if (!this->abcState()) {
     547    if (!solveType)
     548      this->dual(0);
     549    else
     550      this->primal(startUp);
     551  } else {
     552    AbcSimplex * abcModel2=new AbcSimplex(*this);
     553    if (interrupt)
     554      currentAbcModel = abcModel2;
     555    //if (abcSimplex_) {
     556    // move factorization stuff
     557    abcModel2->factorization()->synchronize(this->factorization(),abcModel2);
     558    //}
     559    //abcModel2->startPermanentArrays();
     560    int crashState=abcModel2->abcState()&(256+512+1024);
     561    abcModel2->setAbcState(CLP_ABC_WANTED|crashState);
     562    int ifValuesPass=startUp;
     563    // temp
     564    if (fabs(this->primalTolerance()-1.001e-6)<0.999e-9) {
     565      int type=1;
     566      double diff=this->primalTolerance()-1.001e-6;
     567      printf("Diff %g\n",diff);
     568      if (fabs(diff-1.0e-10)<1.0e-13)
     569        type=2;
     570      else if (fabs(diff-1.0e-11)<1.0e-13)
     571        type=3;
     572#if 0
     573      ClpSimplex * thisModel = static_cast<ClpSimplexOther *> (this)->dualOfModel(1.0,1.0);
     574      if (thisModel) {
     575        printf("Dual of model has %d rows and %d columns\n",
     576               thisModel->numberRows(), thisModel->numberColumns());
     577        thisModel->setOptimizationDirection(1.0);
     578        Idiot info(*thisModel);
     579        info.setStrategy(512 | info.getStrategy());
     580        // Allow for scaling
     581        info.setStrategy(32 | info.getStrategy());
     582        info.setStartingWeight(1.0e3);
     583        info.setReduceIterations(6);
     584        info.crash(50, this->messageHandler(), this->messagesPointer(),false);
     585        // make sure later that primal solution in correct place
     586        // and has correct sign
     587        abcModel2->setupDualValuesPass(thisModel->primalColumnSolution(),
     588                                       thisModel->dualRowSolution(),type);
     589        //thisModel->dual();
     590        delete thisModel;
     591      }
     592#else
     593      if (!solveType) {
     594        this->dual(0);
     595        abcModel2->setupDualValuesPass(this->dualRowSolution(),
     596                                       this->primalColumnSolution(),type);
     597      } else {
     598        ifValuesPass=1;
     599        abcModel2->setStateOfProblem(abcModel2->stateOfProblem() | VALUES_PASS);
     600        Idiot info(*abcModel2);
     601        info.setStrategy(512 | info.getStrategy());
     602        // Allow for scaling
     603        info.setStrategy(32 | info.getStrategy());
     604        info.setStartingWeight(1.0e3);
     605        info.setReduceIterations(6);
     606        info.crash(200, abcModel2->messageHandler(), abcModel2->messagesPointer(),false);
     607        //memcpy(abcModel2->primalColumnSolution(),this->primalColumnSolution(),
     608        //     this->numberColumns()*sizeof(double));
     609      }
     610#endif
     611    }
     612    char line[200];
     613#if ABC_PARALLEL
     614#if ABC_PARALLEL==2
     615#ifndef FAKE_CILK
     616    if (!number_cilk_workers) {
     617      number_cilk_workers=__cilkrts_get_nworkers();
     618      sprintf(line,"%d cilk workers",number_cilk_workers);
     619      handler_->message(CLP_GENERAL, messages_)
     620        << line
     621        << CoinMessageEol;
     622    }
     623#endif
     624#endif
     625    int numberCpu=this->abcState()&15;
     626    if (numberCpu==9) {
     627      numberCpu=1;
     628#if ABC_PARALLEL==2
     629#ifndef FAKE_CILK
     630      if (number_cilk_workers>1)
     631      numberCpu=CoinMin(2*number_cilk_workers,8);
     632#endif
     633#endif
     634    } else if (numberCpu==10) {
     635      // maximum
     636      numberCpu=4;
     637    } else if (numberCpu==10) {
     638      // decide
     639      if (abcModel2->getNumElements()<5000)
     640        numberCpu=1;
     641#if ABC_PARALLEL==2
     642#ifndef FAKE_CILK
     643      else if (number_cilk_workers>1)
     644        numberCpu=CoinMin(2*number_cilk_workers,8);
     645#endif
     646#endif
     647      else
     648        numberCpu=1;
     649    } else {
     650#if ABC_PARALLEL==2
     651#if 0 //ndef FAKE_CILK
     652      char temp[3];
     653      sprintf(temp,"%d",numberCpu);
     654      __cilkrts_set_param("nworkers",temp);
     655#endif
     656#endif
     657    }
     658    abcModel2->setParallelMode(numberCpu-1);
     659#endif
     660    //if (abcState()==3||abcState()==4) {
     661    //abcModel2->setMoreSpecialOptions((65536*2)|abcModel2->moreSpecialOptions());
     662    //}
     663    //if (processTune>0&&processTune<8)
     664    //abcModel2->setMoreSpecialOptions(abcModel2->moreSpecialOptions()|65536*processTune);
     665#if ABC_INSTRUMENT
     666    double startTimeCpu=CoinCpuTime();
     667    double startTimeElapsed=CoinGetTimeOfDay();
     668#if ABC_INSTRUMENT>1
     669    memset(abcPricing,0,sizeof(abcPricing));
     670    memset(abcPricingDense,0,sizeof(abcPricing));
     671    instrument_initialize(abcModel2->numberRows());
     672#endif
     673#endif
     674    if (!solveType)
     675      abcModel2->ClpSimplex::doAbcDual();
     676    else
     677      abcModel2->ClpSimplex::doAbcPrimal(ifValuesPass);
     678#if ABC_INSTRUMENT
     679    double timeCpu=CoinCpuTime()-startTimeCpu;
     680    double timeElapsed=CoinGetTimeOfDay()-startTimeElapsed;
     681    sprintf(line,"Cpu time for %s (%d rows, %d columns %d elements) %g elapsed %g ratio %g - %d iterations",
     682            this->problemName().c_str(),this->numberRows(),this->numberColumns(),
     683            this->getNumElements(),
     684            timeCpu,timeElapsed,timeElapsed ? timeCpu/timeElapsed : 1.0,abcModel2->numberIterations());
     685    handler_->message(CLP_GENERAL, messages_)
     686      << line
     687      << CoinMessageEol;
     688#if ABC_INSTRUMENT>1
     689    {
     690      int n;
     691      n=0;
     692      for (int i=0;i<20;i++)
     693        n+= abcPricing[i];
     694      printf("CCSparse pricing done %d times",n);
     695      int n2=0;
     696      for (int i=0;i<20;i++)
     697        n2+= abcPricingDense[i];
     698      if (n2)
     699        printf(" and dense pricing done %d times\n",n2);
     700      else
     701        printf("\n");
     702      n=0;
     703      printf ("CCS");
     704      for (int i=0;i<19;i++) {
     705        if (abcPricing[i]) {
     706          if (n==5) {
     707            n=0;
     708            printf("\nCCS");
     709          }
     710          n++;
     711          printf("(%d els,%d times) ",i,abcPricing[i]);
     712        }
     713      }
     714      if (abcPricing[19]) {
     715        if (n==5) {
     716          n=0;
     717          printf("\nCCS");
     718        }
     719        n++;
     720        printf("(>=19 els,%d times) ",abcPricing[19]);
     721      }
     722      if (n2) {
     723        printf ("CCD");
     724        for (int i=0;i<19;i++) {
     725          if (abcPricingDense[i]) {
     726            if (n==5) {
     727              n=0;
     728              printf("\nCCD");
     729            }
     730            n++;
     731            int k1=(numberRows_/16)*i;;
     732            int k2=CoinMin(numberRows_,k1+(numberRows_/16)-1);
     733            printf("(%d-%d els,%d times) ",k1,k2,abcPricingDense[i]);
     734          }
     735        }
     736      }
     737      printf("\n");
     738    }
     739    instrument_print();
     740#endif
     741#endif
     742    abcModel2->moveStatusToClp(this);
     743    //ClpModel::stopPermanentArrays();
     744    this->setSpecialOptions(this->specialOptions()&~65536);
     745#if 0
     746    this->writeBasis("a.bas",true);
     747    for (int i=0;i<this->numberRows();i++)
     748      printf("%d %g\n",i,this->dualRowSolution()[i]);
     749    this->dual();
     750    this->writeBasis("b.bas",true);
     751    for (int i=0;i<this->numberRows();i++)
     752      printf("%d %g\n",i,this->dualRowSolution()[i]);
     753#endif
     754    // switch off initialSolve flag
     755    moreSpecialOptions_ &= ~16384;
     756    //this->setNumberIterations(abcModel2->numberIterations()+this->numberIterations());
     757    delete abcModel2;
     758  }
     759}
     760#endif
    360761/** General solve algorithm which can do presolve
    361762    special options (bits)
     
    456857     if (presolve != ClpSolve::presolveOff) {
    457858          bool costedSlacks = false;
     859#ifdef ABC_INHERIT
     860          int numberPasses = 20;
     861#else
    458862          int numberPasses = 5;
     863#endif
    459864          if (presolve == ClpSolve::presolveNumber) {
    460865               numberPasses = options.getPresolvePasses();
     
    501906               presolve = ClpSolve::presolveOff;
    502907          } else {
     908#if 0 //def ABC_INHERIT
     909            {
     910              AbcSimplex * abcModel2=new AbcSimplex(*model2);
     911              delete model2;
     912              model2=abcModel2;
     913              pinfo->setPresolvedModel(model2);
     914            }
     915#else
     916            //ClpModel::stopPermanentArrays();
     917            //setSpecialOptions(specialOptions()&~65536);
     918#endif
    503919              model2->eventHandler()->setSimplex(model2);
    504920              int rcode=model2->eventHandler()->event(ClpEventHandler::presolveSize);
     
    7401156                  // CHECKME This test of specialOptions and the one above
    7411157                  // don't seem compatible.
     1158#ifndef ABC_INHERIT
    7421159                    if ((specialOptions_ & 1024) == 0) {
    7431160                         model2->replaceMatrix(newMatrix);
    7441161                    } else {
    745                          // in integer - just use for sprint/idiot
     1162#endif
     1163                         // in integer (or abc) - just use for sprint/idiot
    7461164                         saveMatrix = NULL;
    7471165                         delete newMatrix;
    748                     }
     1166#ifndef ABC_INHERIT
     1167                    }
     1168#endif
    7491169               } else {
    7501170                    handler_->message(CLP_MATRIX_CHANGE, messages_)
     
    7901210                              doSprint = 0; // switch off sprint
    7911211                    }
    792                     // switch off sprint or idiot if any free variable
     1212                    // switch off idiot or sprint if any free variable
     1213                    // switch off sprint if very few with costs
    7931214                    int iColumn;
    794                     double * columnLower = model2->columnLower();
    795                     double * columnUpper = model2->columnUpper();
     1215                    const double * columnLower = model2->columnLower();
     1216                    const double * columnUpper = model2->columnUpper();
     1217                    const double * objective = model2->objective();
     1218                    int nObj=0;
    7961219                    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
    7971220                         if (columnLower[iColumn] < -1.0e10 && columnUpper[iColumn] > 1.0e10) {
    798                               doSprint = 0;
     1221                              doSprint = 0;
    7991222                              doIdiot = 0;
    8001223                              break;
     1224                         } else if (objective[iColumn]) {
     1225                           nObj++;
    8011226                         }
    8021227                    }
     1228                    if (nObj*10<numberColumns)
     1229                      doSprint=0;
    8031230                    int nPasses = 0;
    8041231                    // look at rhs
     
    8311258                         if (value2 > value1) {
    8321259                              numberNotE++;
    833                               if (value2 > 1.0e31 || value1 < -1.0e31)
    834                                    largestGap = COIN_DBL_MAX;
    835                               else
    836                                    largestGap = value2 - value1;
     1260                              //if (value2 > 1.0e31 || value1 < -1.0e31)
     1261                              //   largestGap = COIN_DBL_MAX;
     1262                              //else
     1263                              //   largestGap = value2 - value1;
    8371264                         }
    8381265                    }
     
    10201447#endif
    10211448          if (doCrash) {
     1449#ifdef ABC_INHERIT
     1450            if (!model2->abcState()) {
     1451#endif
    10221452               switch(doCrash) {
    10231453                    // standard
     
    10351465                    break;
    10361466               }
     1467#ifdef ABC_INHERIT
     1468            } else if (doCrash>=0) {
     1469               model2->setAbcState(model2->abcState()|256*doCrash);
     1470            }
     1471#endif
    10371472          }
    10381473          if (!nPasses) {
     
    10531488                         clpMatrix->releaseSpecialColumnCopy();
    10541489                    } else {
    1055                          model2->dual(0);
     1490#ifndef ABC_INHERIT
     1491                      model2->dual(0);
     1492#else
     1493                      model2->dealWithAbc(0,0,interrupt);
     1494#endif
    10561495                    }
    10571496               } else {
     
    12581697                         info.setStartingWeight(1.0e3);
    12591698                         info.setReduceIterations(6);
    1260                          if (nPasses >= 5000) {
     1699                         //if (nPasses > 200)
     1700                         //info.setFeasibilityTolerance(1.0e-9);
     1701                         //if (nPasses > 1900)
     1702                         //info.setWeightFactor(0.93);
     1703                         if (nPasses > 900) {
     1704                           double reductions=nPasses/6.0;
     1705                           if (nPasses<5000) {
     1706                             reductions /= 12.0;
     1707                           } else {
     1708                             reductions /= 13.0;
     1709                             info.setStartingWeight(1.0e4);
     1710                           }
     1711                           double ratio=1.0/exp10(1.0/reductions);
     1712                           printf("%d passes reduction factor %g\n",nPasses,ratio);
     1713                           info.setWeightFactor(ratio);
     1714                         } else if (nPasses > 500) {
     1715                           info.setWeightFactor(0.7);
     1716                         } else if (nPasses > 200) {
     1717                           info.setWeightFactor(0.5);
     1718                         }
     1719                         if (maximumIterations()<nPasses) {
     1720                           printf("Presuming maximumIterations is just for Idiot\n");
     1721                           nPasses=maximumIterations();
     1722                           setMaximumIterations(COIN_INT_MAX);
     1723                           model2->setMaximumIterations(COIN_INT_MAX);
     1724                         }
     1725                         if (nPasses >= 10000&&nPasses<100000) {
    12611726                              int k = nPasses % 100;
    1262                               nPasses /= 100;
     1727                              nPasses /= 200;
    12631728                              info.setReduceIterations(3);
    12641729                              if (k)
     
    13151780#ifndef LACI_TRY
    13161781                    //if (doIdiot>0)
     1782#ifdef ABC_INHERIT
     1783                    if (!model2->abcState())
     1784#endif
    13171785                    info.setStrategy(512 | info.getStrategy());
    13181786#endif
     
    13271795                              << CoinMessageEol;
    13281796                    timeX = time2;
     1797                    if (nPasses>100000&&nPasses<100500) {
     1798                      // make sure no status left
     1799                      model2->createStatus();
     1800                      // solve
     1801                      if (model2->factorizationFrequency() == 200) {
     1802                        // User did not touch preset
     1803                        model2->defaultFactorizationFrequency();
     1804                      }
     1805                      //int numberRows = model2->numberRows();
     1806                      int numberColumns = model2->numberColumns();
     1807                      // save duals
     1808                      //double * saveDuals = CoinCopyOfArray(model2->dualRowSolution(),numberRows);
     1809                      // for moment this only works on nug etc (i.e. all ==)
     1810                      // needs row objective
     1811                      double * saveObj = CoinCopyOfArray(model2->objective(),numberColumns);
     1812                      double * pi = model2->dualRowSolution();
     1813                      model2->clpMatrix()->transposeTimes(-1.0, pi, model2->objective());
     1814                      // just primal values pass
     1815                      double saveScale = model2->objectiveScale();
     1816                      model2->setObjectiveScale(1.0e-3);
     1817                      model2->primal(2);
     1818                      model2->writeMps("xx.mps");
     1819                      double * solution = model2->primalColumnSolution();
     1820                      double * upper = model2->columnUpper();
     1821                      for (int i=0;i<numberColumns;i++) {
     1822                        if (solution[i]<100.0)
     1823                          upper[i]=1000.0;
     1824                      }
     1825                      model2->setProblemStatus(-1);
     1826                      model2->setObjectiveScale(saveScale);
     1827#ifdef ABC_INHERIT
     1828                      AbcSimplex * abcModel2=new AbcSimplex(*model2);
     1829                      if (interrupt)
     1830                        currentAbcModel = abcModel2;
     1831                      if (abcSimplex_) {
     1832                        // move factorization stuff
     1833                        abcModel2->factorization()->synchronize(model2->factorization(),abcModel2);
     1834                      }
     1835                      abcModel2->startPermanentArrays();
     1836                      abcModel2->setAbcState(CLP_ABC_WANTED);
     1837#if ABC_PARALLEL
     1838                      int parallelMode=1;
     1839                      printf("Parallel mode %d\n",parallelMode);
     1840                      abcModel2->setParallelMode(parallelMode);
     1841#endif
     1842                      //if (processTune>0&&processTune<8)
     1843                      //abcModel2->setMoreSpecialOptions(abcModel2->moreSpecialOptions()|65536*processTune);
     1844                      abcModel2->doAbcDual();
     1845                      abcModel2->moveStatusToClp(model2);
     1846                      //ClpModel::stopPermanentArrays();
     1847                      model2->setSpecialOptions(model2->specialOptions()&~65536);
     1848                      //model2->dual();
     1849                      //model2->setNumberIterations(abcModel2->numberIterations()+model2->numberIterations());
     1850                      delete abcModel2;
     1851#endif
     1852                      memcpy(model2->objective(),saveObj,numberColumns*sizeof(double));
     1853                      //delete [] saveDuals;
     1854                      delete [] saveObj;
     1855                      model2->dual(2);
     1856                    } // end dubious idiot
    13291857               }
    13301858          }
     
    13701898                         clpMatrix->releaseSpecialColumnCopy();
    13711899                    } else {
    1372                          model2->primal(primalStartup);
     1900#ifndef ABC_INHERIT
     1901                        model2->primal(primalStartup);
     1902#else
     1903                        model2->dealWithAbc(1,primalStartup,interrupt);
     1904#endif
    13731905                    }
    13741906               } else {
     1907#ifndef ABC_INHERIT
    13751908                    model2->primal(primalStartup);
     1909#else
     1910                    model2->dealWithAbc(1,primalStartup,interrupt);
     1911#endif
    13761912               }
    13771913          }
     
    16852221               }
    16862222               // Get objective offset
     2223               const double * objective = model2->objective();
    16872224               double offset = 0.0;
    1688                const double * objective = model2->objective();
    1689                for (iColumn = 0; iColumn < numberColumns; iColumn++)
     2225               for (iColumn = 0; iColumn < originalNumberColumns; iColumn++)
    16902226                    offset += fullSolution[iColumn] * objective[iColumn];
     2227#if 0
     2228               // Set artificials to zero if first time close to zero
     2229               for (iColumn = originalNumberColumns; iColumn < numberColumns; iColumn++) {
     2230                 if (fullSolution[iColumn]<primalTolerance_&&objective[iColumn]==penalty) {
     2231                   model2->objective()[iColumn]=2.0*penalty;
     2232                   fullSolution[iColumn]=0.0;
     2233                 }
     2234               }
     2235#endif
    16912236               small.setDblParam(ClpObjOffset, originalOffset - offset);
    16922237               model2->clpMatrix()->times(1.0, fullSolution, sumFixed);
     
    17172262                    } else {
    17182263#if 1
    1719                          small.primal(1);
     2264#ifdef ABC_INHERIT
     2265                      //small.writeMps("try.mps");
     2266                      if (iPass)
     2267                         small.dealWithAbc(1,1);
     2268                       else
     2269                         small.dealWithAbc(0,0);
     2270#else
     2271                      if (iPass)
     2272                         small.primal(1);
     2273                      else
     2274                         small.dual(0);
     2275#endif
    17202276#else
    17212277                         int numberColumns = small.numberColumns();
     
    17292285                                               nBound, false, false);
    17302286                         if (small2) {
    1731                               small2->primal(1);
     2287#ifdef ABC_INHERIT
     2288                              small2->dealWithAbc(1,1);
     2289#else
     2290                              small.primal(1);
     2291#endif
    17322292                              if (small2->problemStatus() == 0) {
    17332293                                   small.setProblemStatus(0);
    17342294                                   ((ClpSimplexOther *) (&small))->afterCrunch(*small2, whichRow, whichColumn, nBound);
    17352295                              } else {
    1736                                    small2->primal(1);
     2296#ifdef ABC_INHERIT
     2297                                   small2->dealWithAbc(1,1);
     2298#else
     2299                                   small.primal(1);
     2300#endif
    17372301                                   if (small2->problemStatus())
    17382302                                        small.primal(1);
     
    18692433               delete [] saveUpper;
    18702434          }
     2435#ifdef ABC_INHERIT
     2436          model2->dealWithAbc(1,1);
     2437#else
    18712438          model2->primal(1);
     2439#endif
    18722440          model2->setPerturbation(savePerturbation);
    18732441          if (model2 != originalModel2) {
     
    18982466          if (interrupt)
    18992467               currentModel2 = &barrier;
     2468          if (barrier.numberRows()+barrier.numberColumns()>10000)
     2469            barrier.setMaximumBarrierIterations(1000);
    19002470          int barrierOptions = options.getSpecialOption(4);
    19012471          int aggressiveGamma = 0;
     
    19412511          default:
    19422512               if (!doKKT) {
    1943                     ClpCholeskyBase * cholesky = new ClpCholeskyBase();
     2513                    ClpCholeskyBase * cholesky = new ClpCholeskyBase(options.getExtraInfo(1));
    19442514                    cholesky->setIntegerParameter(0, speed);
    19452515                    barrier.setCholesky(cholesky);
     
    22612831                         model2->defaultFactorizationFrequency();
    22622832                    }
    2263 #if 1
     2833#if 1 //ndef ABC_INHERIT //#if 1
    22642834                    // throw some into basis
    22652835                    if(!forceFixing) {
     
    23202890                              int numberRows = model2->numberRows();
    23212891                              int numberColumns = model2->numberColumns();
     2892#ifdef ABC_INHERIT
     2893                              model2->checkSolution(0);
     2894                              printf("%d primal infeasibilities summing to %g\n",
     2895                                     model2->numberPrimalInfeasibilities(),
     2896                                     model2->sumPrimalInfeasibilities());
     2897                              model2->dealWithAbc(1,1);
     2898                         }
     2899                    }
     2900#else
    23222901                              // just primal values pass
    23232902                              double saveScale = model2->objectiveScale();
     
    24002979                    model2->setObjectiveScale(saveScale);
    24012980                    model2->primal(1);
     2981#endif
    24022982#else
    24032983                    // just primal
    2404                     model2->primal(1);
     2984#ifdef ABC_INHERIT
     2985                    model2->checkSolution(0);
     2986                    printf("%d primal infeasibilities summing to %g\n",
     2987                           model2->numberPrimalInfeasibilities(),
     2988                           model2->sumPrimalInfeasibilities());
     2989                    model2->dealWithAbc(1,1);
     2990#else
     2991                    model2->primal(1);
     2992#endif
     2993                    //model2->primal(1);
    24052994#endif
    24062995               } else if (barrierStatus == 4) {
     
    24893078                    << CoinMessageEol;
    24903079          timeX = time2;
    2491           if (!presolveToFile)
     3080          if (!presolveToFile) {
     3081#if 1 //ndef ABC_INHERIT
    24923082               delete model2;
     3083#else
     3084               if (model2->abcSimplex())
     3085                 delete model2->abcSimplex();
     3086               else
     3087                 delete model2;
     3088#endif
     3089          }
    24933090          if (interrupt)
    24943091               currentModel = this;
     
    25123109                         dual();
    25133110                    } else {
    2514                          setPerturbation(100);
    2515                          primal(1);
     3111                      if (numberRows_<10000)
     3112                        setPerturbation(100); // probably better to perturb after n its
     3113                      else if (savePerturbation<100)
     3114                        setPerturbation(51); // probably better to perturb after n its
     3115#ifndef ABC_INHERIT
     3116                        primal(1);
     3117#else
     3118                        dealWithAbc(1,1,interrupt);
     3119#endif
    25163120                    }
    25173121               } else {
     
    25293133               timeX = time2;
    25303134          } else if (rcode >= 0) {
    2531               primal(1);
     3135#ifdef ABC_INHERIT
     3136            dealWithAbc(1,1,true);
     3137#else
     3138            primal(1);
     3139#endif
    25323140          } else {
    25333141               secondaryStatus_ = finalSecondaryStatus;
     
    29623570     if (!model_)
    29633571          return -1;
    2964      double objective = model_->rawObjectiveValue();
    2965      if (model_->algorithm() < 0)
     3572     double objective;
     3573     if (model_->algorithm() < 0) {
     3574       objective = model_->rawObjectiveValue();
    29663575          objective -= model_->bestPossibleImprovement();
     3576     } else {
     3577       objective = model_->rawObjectiveValue();
     3578     }
    29673579     double infeasibility;
    29683580     double realInfeasibility = 0.0;
  • trunk/Clp/src/ClpSolve.hpp

    r1665 r1878  
    385385     //@}
    386386};
     387#if CLP_HAS_ABC
     388#include "AbcCommon.hpp"
     389/// For saving extra information to see if looping.
     390class AbcSimplexProgress : public ClpSimplexProgress {
     391
     392public:
     393
     394
     395     /**@name Constructors and destructor and copy */
     396     //@{
     397     /// Default constructor
     398     AbcSimplexProgress (  );
     399
     400     /// Constructor from model
     401     AbcSimplexProgress ( ClpSimplex * model );
     402
     403     /// Copy constructor.
     404     AbcSimplexProgress(const AbcSimplexProgress &);
     405
     406     /// Assignment operator. This copies the data
     407     AbcSimplexProgress & operator=(const AbcSimplexProgress & rhs);
     408     /// Destructor
     409     ~AbcSimplexProgress (  );
     410
     411     //@}
     412
     413     /**@name Check progress */
     414     //@{
     415     /** Returns -1 if okay, -n+1 (n number of times bad) if bad but action taken,
     416         >=0 if give up and use as problem status
     417     */
     418     int looping (  );
     419
     420     //@}
     421     /**@name Data  */
     422     //@}
     423};
    387424#endif
     425#endif
  • trunk/Clp/src/IdiSolve.cpp

    r1723 r1878  
    130130IdiotResult
    131131Idiot::IdiSolve(
    132      int nrows, int ncols, double * rowsol , double * colsol,
    133      double * pi, double * djs, const double * origcost , double * rowlower,
    134      double * rowupper, const double * lower,
    135      const double * upper, const double * elemnt,
     132     int nrows, int ncols, double * COIN_RESTRICT rowsol , double * COIN_RESTRICT colsol,
     133     double * COIN_RESTRICT pi, double * COIN_RESTRICT djs, const double * COIN_RESTRICT origcost , double * COIN_RESTRICT rowlower,
     134     double * COIN_RESTRICT rowupper, const double * COIN_RESTRICT lower,
     135     const double * COIN_RESTRICT upper, const double * COIN_RESTRICT elemnt,
    136136     const int * row, const CoinBigIndex * columnStart,
    137      const int * length, double * lambda,
     137     const int * length, double * COIN_RESTRICT lambda,
    138138     int maxIts, double mu, double drop,
    139139     double maxmin, double offset,
     
    150150     int extraBlock = 0;
    151151     int * rowExtra = NULL;
    152      double * solExtra = NULL;
    153      double * elemExtra = NULL;
    154      double * upperExtra = NULL;
    155      double * costExtra = NULL;
    156      double * useCostExtra = NULL;
    157      double * saveExtra = NULL;
    158      double * cost = NULL;
     152     double * COIN_RESTRICT solExtra = NULL;
     153     double * COIN_RESTRICT elemExtra = NULL;
     154     double * COIN_RESTRICT upperExtra = NULL;
     155     double * COIN_RESTRICT costExtra = NULL;
     156     double * COIN_RESTRICT useCostExtra = NULL;
     157     double * COIN_RESTRICT saveExtra = NULL;
     158     double * COIN_RESTRICT cost = NULL;
    159159     double saveValue = 1.0e30;
    160160     double saveOffset = offset;
     
    167167#endif
    168168     int nflagged;
    169      double * thetaX;
    170      double * djX;
    171      double * bX;
    172      double * vX;
     169     double * COIN_RESTRICT thetaX;
     170     double * COIN_RESTRICT djX;
     171     double * COIN_RESTRICT bX;
     172     double * COIN_RESTRICT vX;
    173173     double ** aX;
    174174     double **aworkX;
    175175     double ** allsum;
    176      double * saveSol = 0;
    177      const double * useCost = cost;
     176     double * COIN_RESTRICT saveSol = 0;
     177     const double * COIN_RESTRICT useCost = cost;
    178178     double bestSol = 1.0e60;
    179179     double weight = 0.5 / mu;
     
    197197          obj[i] = 1.0e70;
    198198     }
    199      //#define TWO_GOES
    200 #ifdef TWO_GOES
    201      double * pi2 = new double [nrows];
    202      double * rowsol2 = new double [nrows];
     199     //#define FOUR_GOES 2
     200#ifdef FOUR_GOES
     201     double * COIN_RESTRICT pi2 = new double [3*nrows];
     202     double * COIN_RESTRICT rowsol2 = new double [3*nrows];
     203     double * COIN_RESTRICT piX[4];
     204     double * COIN_RESTRICT rowsolX[4];
     205     int startsX[2][5];
     206     int nChangeX[4];
     207     double maxDjX[4];
     208     double objvalueX[4];
     209     int nflaggedX[4];
     210     piX[0]=pi;
     211     piX[1]=pi2;
     212     piX[2]=pi2+nrows;
     213     piX[3]=piX[2]+nrows;
     214     rowsolX[0]=rowsol;
     215     rowsolX[1]=rowsol2;
     216     rowsolX[2]=rowsol2+nrows;
     217     rowsolX[3]=rowsolX[2]+nrows;
    203218#endif
    204219     allsum = new double * [nsolve];
     
    339354     iter = 0;
    340355     for (; iter < maxIts; iter++) {
    341           double sum1 = 0.0, sum2 = 0.0, djtot;
     356          double sum1 = 0.0, sum2 = 0.0;
    342357          double lastObj = 1.0e70;
    343358          int good = 0, doScale = 0;
     
    346361               ii = ii * EVERY;
    347362               if (iter > ii - HISTORY * 2 && (iter & 1) == 0) {
    348                     double * x = history[HISTORY-1];
     363                    double * COIN_RESTRICT x = history[HISTORY-1];
    349364                    for (i = HISTORY - 1; i > 0; i--) {
    350365                         history[i] = history[i-1];
     
    357372          if ((iter % SAVEHISTORY) == 0 || doFull) {
    358373               if ((strategy & 16) == 0) {
    359                     double * x = history[HISTORY-1];
     374                    double * COIN_RESTRICT x = history[HISTORY-1];
    360375                    for (i = HISTORY - 1; i > 0; i--) {
    361376                         history[i] = history[i-1];
     
    380395                    start[1] = 0;
    381396                    stop[1] = kcol;
     397#ifdef FOUR_GOES
     398                    for (int itry=0;itry<2;itry++) {
     399                      int chunk=(stop[itry]-start[itry]+FOUR_GOES-1)/FOUR_GOES;
     400                      startsX[itry][0]=start[itry];
     401                      for (int i=1;i<5;i++)
     402                        startsX[itry][i]=CoinMin(stop[itry],startsX[itry][i-1]+chunk);
     403                    }
     404#endif
    382405               } else {
    383406                    start[0] = kcol;
     
    385408                    start[1] = ncols - 1;
    386409                    stop[1] = kcol;
     410#ifdef FOUR_GOES
     411                    for (int itry=0;itry<2;itry++) {
     412                      int chunk=(start[itry]-stop[itry]+FOUR_GOES-1)/FOUR_GOES;
     413                      startsX[itry][0]=start[itry];
     414                      for (int i=1;i<5;i++)
     415                        startsX[itry][i]=CoinMax(stop[itry],startsX[itry][i-1]-chunk);
     416                    }
     417#endif
    387418               }
    388419               int itry = 0;
    389420               /*if ((strategy&16)==0) {
    390                 double * x=history[HISTORY-1];
     421                double * COIN_RESTRICT x=history[HISTORY-1];
    391422                for (i=HISTORY-1;i>0;i--) {
    392423                history[i]=history[i-1];
     
    407438                         objvalue = 0.0;
    408439                         memset(pi, 0, nrows * sizeof(double));
    409                          djtot = 0.0;
    410440                         {
    411                               double * theta = thetaX;
    412                               double * dj = djX;
    413                               double * b = bX;
     441                              double * COIN_RESTRICT theta = thetaX;
     442                              double * COIN_RESTRICT dj = djX;
     443                              double * COIN_RESTRICT b = bX;
    414444                              double ** a = aX;
    415445                              double ** awork = aworkX;
    416                               double * v = vX;
     446                              double * COIN_RESTRICT v = vX;
    417447                              double c;
    418448#ifdef FIT
     
    883913          maxDj = 0.0;
    884914          // go through forwards or backwards and starting at odd places
    885           int itry;
    886 #ifdef TWO_GOES
    887           memcpy(pi2, pi, nrows * sizeof(double));
    888           memcpy(rowsol2, rowsol, nrows * sizeof(double));
    889 #endif
    890           for (itry = 0; itry < 2; itry++) {
    891                int icol = start[itry];
     915#ifdef FOUR_GOES
     916          for (int i=1;i<FOUR_GOES;i++) {
     917            cilk_spawn memcpy(piX[i], pi, nrows * sizeof(double));
     918            cilk_spawn memcpy(rowsolX[i], rowsol, nrows * sizeof(double));
     919          }
     920          for (int i=0;i<FOUR_GOES;i++) {
     921            nChangeX[i]=0;
     922            maxDjX[i]=0.0;
     923            objvalueX[i]=0.0;
     924            nflaggedX[i]=0;
     925          }
     926          cilk_sync;
     927#endif
     928          //printf("PASS\n");
     929#ifdef FOUR_GOES
     930               cilk_for (int iPar = 0; iPar < FOUR_GOES; iPar++) {
     931                    double * COIN_RESTRICT pi=piX[iPar];
     932                    double * COIN_RESTRICT rowsol = rowsolX[iPar];
     933          for (int itry = 0; itry < 2; itry++) {
     934                    int istop;
     935                    int istart;
     936#if 0
     937                    int chunk = (start[itry]+stop[itry]+FOUR_GOES-1)/FOUR_GOES;
     938                    if (iPar == 0) {
     939                      istart=start[itry];
     940                      istop=(start[itry]+stop[itry])>>1;
     941                    } else {
     942                      istart=(start[itry]+stop[itry])>>1;
     943                      istop = stop[itry];
     944                    }
     945#endif
     946#if 0
     947                    printf("istart %d istop %d direction %d array %d %d new %d %d\n",
     948                           istart,istop,direction,start[itry],stop[itry],
     949                           startsX[itry][iPar],startsX[itry][iPar+1]);
     950#endif
     951                    istart=startsX[itry][iPar];
     952                    istop=startsX[itry][iPar+1];
     953#else
     954          for (int itry = 0; itry < 2; itry++) {
     955               int istart = start[itry];
    892956               int istop = stop[itry];
    893 #ifdef TWO_GOES
    894                for (int iPar = 0; iPar < 2; iPar++) {
    895                     double * temp = pi;
    896                     pi = pi2;
    897                     pi2 = temp;
    898                     temp = rowsol;
    899                     rowsol = rowsol2;
    900                     rowsol2 = temp;
    901                     if (iPar == 0) {
    902                          istop = (icol + istop) >> 1;
    903                     } else {
    904                          icol = istop;
    905                          istop = stop[itry];
    906                     }
    907 #endif
    908                     for (; icol != istop; icol += direction) {
     957#endif
     958                    for (int icol=istart; icol != istop; icol += direction) {
    909959                         if (!statusWork[icol]) {
    910960                              CoinBigIndex j;
     
    935985                              if (djval > djTol) {
    936986                                   if (djval2 < -1.0e-4) {
    937                                         double valuep, thetap;
     987#ifndef FOUR_GOES
    938988                                        nChange++;
    939989                                        if (djval > maxDj) maxDj = djval;
     990#else
     991                                        nChangeX[iPar]++;
     992                                        if (djval > maxDjX[iPar]) maxDjX[iPar] = djval;
     993#endif
    940994                                        /*if (djval>3.55e6) {
    941995                                                        printf("big\n");
     
    9671021                                        /* solve */
    9681022                                        theta = -b / a;
     1023#ifndef FOUR_GOES
    9691024                                        if ((strategy & 4) != 0) {
     1025                                          double valuep, thetap;
    9701026                                             value2 = a * theta * theta + 2.0 * b * theta;
    9711027                                             thetap = 2.0 * theta;
     
    9781034                                             }
    9791035                                        }
     1036#endif
    9801037                                        if (theta > 0.0) {
    9811038                                             if (theta < upper[icol] - colsol[icol]) {
     
    9921049                                        }
    9931050                                        colsol[icol] += value2;
     1051#ifndef FOUR_GOES
    9941052                                        objvalue += cost[icol] * value2;
     1053#else
     1054                                        objvalueX[iPar] += cost[icol] * value2;
     1055#endif
    9951056                                        if (elemnt) {
    9961057                                             for (j = columnStart[icol]; j < columnStart[icol] + length[icol]; j++) {
     
    10141075                                        if (djval > djFlag) {
    10151076                                             statusWork[icol] = 1;
     1077#ifndef FOUR_GOES
    10161078                                             nflagged++;
    1017                                         }
    1018                                    }
    1019                               }
    1020                          }
    1021                     }
    1022 #ifdef TWO_GOES
    1023                }
    1024 #endif
    1025           }
    1026 #ifdef TWO_GOES
    1027           for (i = 0; i < nrows; i++) {
    1028                rowsol[i] = 0.5 * (rowsol[i] + rowsol2[i]);
    1029                pi[i] = 0.5 * (pi[i] + pi2[i]);
     1079#else
     1080                                             nflaggedX[iPar]++;
     1081#endif
     1082                                        }
     1083                                   }
     1084                              }
     1085                         }
     1086                    }
     1087#ifdef FOUR_GOES
     1088               }
     1089#endif
     1090          }
     1091#ifdef FOUR_GOES
     1092          for (int i=0;i<FOUR_GOES;i++) {
     1093            nChange += nChangeX[i];
     1094            maxDj = CoinMax(maxDj,maxDjX[i]);
     1095            objvalue += objvalueX[i];
     1096            nflagged += nflaggedX[i];
     1097          }
     1098          cilk_for (int i = 0; i < nrows; i++) {
     1099#if FOUR_GOES==2
     1100            rowsol[i] = 0.5 * (rowsolX[0][i] + rowsolX[1][i]);
     1101            pi[i] = 0.5 * (piX[0][i] + piX[1][i]);
     1102#elif FOUR_GOES==3
     1103            pi[i] = 0.33333333333333 * (piX[0][i] + piX[1][i]+piX[2][i]);
     1104            rowsol[i] = 0.3333333333333 * (rowsolX[0][i] + rowsolX[1][i]+rowsolX[2][i]);
     1105#else
     1106            pi[i] = 0.25 * (piX[0][i] + piX[1][i]+piX[2][i]+piX[3][i]);
     1107            rowsol[i] = 0.25 * (rowsolX[0][i] + rowsolX[1][i]+rowsolX[2][i]+rowsolX[3][i]);
     1108#endif
    10301109          }
    10311110#endif
    10321111          if (extraBlock) {
    1033                for (i = 0; i < extraBlock; i++) {
     1112               for (int i = 0; i < extraBlock; i++) {
    10341113                    double value = solExtra[i];
    10351114                    double djval = costExtra[i];
     
    10751154          }
    10761155          if ((iter % 10) == 2) {
    1077                for (i = DJTEST - 1; i > 0; i--) {
     1156               for (int i = DJTEST - 1; i > 0; i--) {
    10781157                    djSave[i] = djSave[i-1];
    10791158               }
     
    10811160               largestDj = CoinMax(largestDj, maxDj);
    10821161               smallestDj = CoinMin(smallestDj, maxDj);
    1083                for (i = DJTEST - 1; i > 0; i--) {
     1162               for (int i = DJTEST - 1; i > 0; i--) {
    10841163                    maxDj += djSave[i];
    10851164               }
     
    11351214     delete [] allsum;
    11361215     delete [] cost;
    1137 #ifdef TWO_GOES
     1216#ifdef FOUR_GOES
    11381217     delete [] pi2 ;
    11391218     delete [] rowsol2 ;
  • trunk/Clp/src/Idiot.cpp

    r1723 r1878  
    1515#include "CoinMessageHandler.hpp"
    1616#include "CoinHelperFunctions.hpp"
     17#if CLP_HAS_ABC
     18#include "AbcCommon.hpp"
     19#endif
    1720// Redefine stuff for Clp
    1821#ifndef OSI_IDIOT
     
    210213                    double rowValue = rowsol[i];
    211214                    CoinBigIndex j = columnStart[iCol];
    212                     rowSave += (colsol[iCol] - lower[iCol]) * element[j];
    213                     colsol[iCol] = lower[iCol];
    214                     assert (lower[iCol] > -1.0e20);
     215                    double lowerValue = CoinMax(CoinMin(colsol[iCol], 0.0) - 1000.0, lower[iCol]);
     216                    rowSave += (colsol[iCol] - lowerValue) * element[j];
     217                    colsol[iCol] = lowerValue;
    215218                    while (nextSlack[iCol] >= 0) {
    216219                         iCol = nextSlack[iCol];
    217220                         j = columnStart[iCol];
    218                          rowSave += (colsol[iCol] - lower[iCol]) * element[j];
    219                          colsol[iCol] = lower[iCol];
     221                         double lowerValue = CoinMax(CoinMin(colsol[iCol], 0.0) - 1000.0, lower[iCol]);
     222                         rowSave += (colsol[iCol] - lowerValue) * element[j];
     223                         colsol[iCol] = lowerValue;
    220224                    }
    221225                    iCol = negSlack[i];
     
    830834#endif
    831835          }
    832           if (iteration > 50 && n == numberAway && result.infeas < 1.0e-4) {
     836          if (iteration > 50 && n == numberAway ) {
     837            if((result.infeas < 1.0e-4 && majorIterations_<200)||result.infeas<1.0e-8) {
    833838#ifdef CLP_INVESTIGATE
    834839               printf("infeas small %g\n", result.infeas);
    835840#endif
    836841               break; // not much happening
     842            }
    837843          }
    838844          if (lightWeight_ == 1 && iteration > 10 && result.infeas > 1.0 && maxIts != 7) {
     
    874880          }
    875881          bestInfeas = CoinMin(bestInfeas, result.infeas);
     882          if (majorIterations_>100&&majorIterations_<200) {
     883            if (iteration==majorIterations_-100) {
     884              // redo
     885              double muX=mu*10.0;
     886              bestInfeas=1.0e3;
     887              mu=muX;
     888              nTry=0;
     889            }
     890          }
    876891          if (iteration) {
    877892               /* this code is in to force it to terminate sometime */
     
    17021717               if (!wantVector) {
    17031718                    //#define TWO_GOES
     1719#ifdef ABC_INHERIT
     1720#ifndef TWO_GOES
     1721                    model_->dealWithAbc(1,justValuesPass ? 2 : 1);
     1722#else
     1723                    model_->dealWithAbc(1,1 + 11);
     1724#endif
     1725#else
    17041726#ifndef TWO_GOES
    17051727                    model_->primal(justValuesPass ? 2 : 1);
    17061728#else
    17071729                    model_->primal(1 + 11);
     1730#endif
    17081731#endif
    17091732               } else {
     
    17121735                    assert (clpMatrix);
    17131736                    clpMatrix->makeSpecialColumnCopy();
     1737#ifdef ABC_INHERIT
     1738                    model_->dealWithAbc(1,1);
     1739#else
    17141740                    model_->primal(1);
     1741#endif
    17151742                    clpMatrix->releaseSpecialColumnCopy();
    17161743               }
    17171744               if (presolve) {
     1745                 model_->primal();
    17181746                    pinfo.postsolve(true);
    17191747                    delete model_;
     
    17281756               saveModel = NULL;
    17291757               if (justValuesPass)
     1758#ifdef ABC_INHERIT
     1759                    model_->dealWithAbc(1,2);
     1760#else
    17301761                    model_->primal(2);
     1762#endif
    17311763          }
    17321764          if (allowInfeasible) {
     
    17801812               }
    17811813               if (!wantVector) {
     1814#ifdef ABC_INHERIT
     1815                    model_->dealWithAbc(1,1);
     1816#else
    17821817                    model_->primal(1);
     1818#endif
    17831819               } else {
    17841820                    ClpMatrixBase * matrix = model_->clpMatrix();
     
    17861822                    assert (clpMatrix);
    17871823                    clpMatrix->makeSpecialColumnCopy();
     1824#ifdef ABC_INHERIT
     1825                    model_->dealWithAbc(1,1);
     1826#else
    17881827                    model_->primal(1);
     1828#endif
    17891829                    clpMatrix->releaseSpecialColumnCopy();
    17901830               }
     
    18201860               }
    18211861               if (!wantVector) {
     1862#ifdef ABC_INHERIT
     1863                    model_->dealWithAbc(1,1);
     1864#else
    18221865                    model_->primal(1);
     1866#endif
    18231867               } else {
    18241868                    ClpMatrixBase * matrix = model_->clpMatrix();
     
    18261870                    assert (clpMatrix);
    18271871                    clpMatrix->makeSpecialColumnCopy();
     1872#ifdef ABC_INHERIT
     1873                    model_->dealWithAbc(1,1);
     1874#else
    18281875                    model_->primal(1);
     1876#endif
    18291877                    clpMatrix->releaseSpecialColumnCopy();
    18301878               }
  • trunk/Clp/src/Makefile.am

    r1763 r1878  
    6969        ClpSimplexPrimal.cpp ClpSimplexPrimal.hpp \
    7070        ClpSolve.cpp ClpSolve.hpp \
     71        AbcWarmStart.cpp AbcWarmStart.hpp \
     72        CoinAbcBaseFactorization.hpp \
     73        CoinAbcCommonFactorization.hpp \
     74        CoinAbcCommon.hpp \
     75        AbcDualRowDantzig.hpp \
     76        AbcDualRowPivot.hpp \
     77        AbcDualRowSteepest.hpp \
     78        CoinAbcFactorization.hpp \
     79        CoinAbcHelperFunctions.hpp \
     80        AbcMatrix.hpp \
     81        AbcCommon.hpp \
     82        AbcSimplexDual.hpp \
     83        AbcSimplexPrimal.hpp \
     84        AbcSimplex.hpp \
     85        CoinAbcBaseFactorization1.cpp \
     86        CoinAbcBaseFactorization2.cpp \
     87        CoinAbcBaseFactorization3.cpp \
     88        CoinAbcBaseFactorization4.cpp \
     89        CoinAbcBaseFactorization5.cpp \
     90        CoinAbcDenseFactorization.cpp \
     91        CoinAbcDenseFactorization.hpp \
     92        AbcDualRowDantzig.cpp \
     93        AbcDualRowPivot.cpp \
     94        AbcDualRowSteepest.cpp \
     95        CoinAbcFactorization1.cpp \
     96        CoinAbcFactorization2.cpp \
     97        CoinAbcFactorization3.cpp \
     98        CoinAbcFactorization4.cpp \
     99        CoinAbcFactorization5.cpp \
     100        CoinAbcHelperFunctions.cpp \
     101        AbcMatrix.cpp \
     102        AbcSimplex.cpp \
     103        AbcSimplexDual.cpp \
     104        AbcSimplexPrimal.cpp \
     105        AbcSimplexParallel.cpp \
     106        AbcPrimalColumnDantzig.cpp AbcPrimalColumnDantzig.hpp \
     107        AbcPrimalColumnPivot.cpp AbcPrimalColumnPivot.hpp \
     108        AbcPrimalColumnSteepest.cpp AbcPrimalColumnSteepest.hpp \
     109        CoinAbcSmallFactorization1.cpp \
     110        CoinAbcSmallFactorization2.cpp \
     111        CoinAbcSmallFactorization3.cpp \
     112        CoinAbcSmallFactorization4.cpp \
     113        CoinAbcSmallFactorization5.cpp \
     114        CoinAbcOrderedFactorization1.cpp \
     115        CoinAbcOrderedFactorization2.cpp \
     116        CoinAbcOrderedFactorization3.cpp \
     117        CoinAbcOrderedFactorization4.cpp \
     118        CoinAbcOrderedFactorization5.cpp \
     119        AbcSimplexFactorization.cpp \
     120        AbcSimplexFactorization.hpp \
     121        AbcNonLinearCost.cpp AbcNonLinearCost.hpp \
    71122        Idiot.cpp Idiot.hpp \
    72123        IdiSolve.cpp
     
    185236        ClpSolve.hpp \
    186237        CbcOrClpParam.hpp \
     238        AbcSimplex.hpp CoinAbcCommon.hpp AbcCommon.hpp \
    187239        Idiot.hpp
    188240
  • trunk/Clp/src/Makefile.in

    r1763 r1878  
    118118        ClpSimplexNonlinear.hpp ClpSimplexOther.cpp \
    119119        ClpSimplexOther.hpp ClpSimplexPrimal.cpp ClpSimplexPrimal.hpp \
    120         ClpSolve.cpp ClpSolve.hpp Idiot.cpp Idiot.hpp IdiSolve.cpp \
     120        ClpSolve.cpp ClpSolve.hpp AbcWarmStart.cpp AbcWarmStart.hpp \
     121        CoinAbcBaseFactorization.hpp CoinAbcCommonFactorization.hpp \
     122        CoinAbcCommon.hpp AbcDualRowDantzig.hpp AbcDualRowPivot.hpp \
     123        AbcDualRowSteepest.hpp CoinAbcFactorization.hpp \
     124        CoinAbcHelperFunctions.hpp AbcMatrix.hpp AbcCommon.hpp \
     125        AbcSimplexDual.hpp AbcSimplexPrimal.hpp AbcSimplex.hpp \
     126        CoinAbcBaseFactorization1.cpp CoinAbcBaseFactorization2.cpp \
     127        CoinAbcBaseFactorization3.cpp CoinAbcBaseFactorization4.cpp \
     128        CoinAbcBaseFactorization5.cpp CoinAbcDenseFactorization.cpp \
     129        CoinAbcDenseFactorization.hpp AbcDualRowDantzig.cpp \
     130        AbcDualRowPivot.cpp AbcDualRowSteepest.cpp \
     131        CoinAbcFactorization1.cpp CoinAbcFactorization2.cpp \
     132        CoinAbcFactorization3.cpp CoinAbcFactorization4.cpp \
     133        CoinAbcFactorization5.cpp CoinAbcHelperFunctions.cpp \
     134        AbcMatrix.cpp AbcSimplex.cpp AbcSimplexDual.cpp \
     135        AbcSimplexPrimal.cpp AbcSimplexParallel.cpp \
     136        AbcPrimalColumnDantzig.cpp AbcPrimalColumnDantzig.hpp \
     137        AbcPrimalColumnPivot.cpp AbcPrimalColumnPivot.hpp \
     138        AbcPrimalColumnSteepest.cpp AbcPrimalColumnSteepest.hpp \
     139        CoinAbcSmallFactorization1.cpp CoinAbcSmallFactorization2.cpp \
     140        CoinAbcSmallFactorization3.cpp CoinAbcSmallFactorization4.cpp \
     141        CoinAbcSmallFactorization5.cpp \
     142        CoinAbcOrderedFactorization1.cpp \
     143        CoinAbcOrderedFactorization2.cpp \
     144        CoinAbcOrderedFactorization3.cpp \
     145        CoinAbcOrderedFactorization4.cpp \
     146        CoinAbcOrderedFactorization5.cpp AbcSimplexFactorization.cpp \
     147        AbcSimplexFactorization.hpp AbcNonLinearCost.cpp \
     148        AbcNonLinearCost.hpp Idiot.cpp Idiot.hpp IdiSolve.cpp \
    121149        ClpCholeskyUfl.cpp ClpCholeskyUfl.hpp ClpCholeskyMumps.cpp \
    122150        ClpCholeskyMumps.hpp ClpCholeskyWssmp.cpp ClpCholeskyWssmp.hpp \
     
    145173        ClpQuadraticObjective.lo ClpSimplex.lo ClpSimplexDual.lo \
    146174        ClpSimplexNonlinear.lo ClpSimplexOther.lo ClpSimplexPrimal.lo \
    147         ClpSolve.lo Idiot.lo IdiSolve.lo $(am__objects_1) \
     175        ClpSolve.lo AbcWarmStart.lo CoinAbcBaseFactorization1.lo \
     176        CoinAbcBaseFactorization2.lo CoinAbcBaseFactorization3.lo \
     177        CoinAbcBaseFactorization4.lo CoinAbcBaseFactorization5.lo \
     178        CoinAbcDenseFactorization.lo AbcDualRowDantzig.lo \
     179        AbcDualRowPivot.lo AbcDualRowSteepest.lo \
     180        CoinAbcFactorization1.lo CoinAbcFactorization2.lo \
     181        CoinAbcFactorization3.lo CoinAbcFactorization4.lo \
     182        CoinAbcFactorization5.lo CoinAbcHelperFunctions.lo \
     183        AbcMatrix.lo AbcSimplex.lo AbcSimplexDual.lo \
     184        AbcSimplexPrimal.lo AbcSimplexParallel.lo \
     185        AbcPrimalColumnDantzig.lo AbcPrimalColumnPivot.lo \
     186        AbcPrimalColumnSteepest.lo CoinAbcSmallFactorization1.lo \
     187        CoinAbcSmallFactorization2.lo CoinAbcSmallFactorization3.lo \
     188        CoinAbcSmallFactorization4.lo CoinAbcSmallFactorization5.lo \
     189        CoinAbcOrderedFactorization1.lo \
     190        CoinAbcOrderedFactorization2.lo \
     191        CoinAbcOrderedFactorization3.lo \
     192        CoinAbcOrderedFactorization4.lo \
     193        CoinAbcOrderedFactorization5.lo AbcSimplexFactorization.lo \
     194        AbcNonLinearCost.lo Idiot.lo IdiSolve.lo $(am__objects_1) \
    148195        $(am__objects_2) $(am__objects_3) $(am__objects_4) \
    149196        $(am__objects_5)
     
    197244        ClpQuadraticObjective.hpp ClpSimplex.hpp \
    198245        ClpSimplexNonlinear.hpp ClpSimplexOther.hpp ClpSimplexDual.hpp \
    199         ClpSimplexPrimal.hpp ClpSolve.hpp CbcOrClpParam.hpp Idiot.hpp \
     246        ClpSimplexPrimal.hpp ClpSolve.hpp CbcOrClpParam.hpp \
     247        AbcSimplex.hpp CoinAbcCommon.hpp AbcCommon.hpp Idiot.hpp \
    200248        ClpCholeskyUfl.hpp ClpCholeskyMumps.hpp ClpCholeskyWssmp.hpp \
    201249        ClpCholeskyWssmpKKT.hpp CbcOrClpParam.cpp
     
    501549        ClpSimplexNonlinear.hpp ClpSimplexOther.cpp \
    502550        ClpSimplexOther.hpp ClpSimplexPrimal.cpp ClpSimplexPrimal.hpp \
    503         ClpSolve.cpp ClpSolve.hpp Idiot.cpp Idiot.hpp IdiSolve.cpp \
     551        ClpSolve.cpp ClpSolve.hpp AbcWarmStart.cpp AbcWarmStart.hpp \
     552        CoinAbcBaseFactorization.hpp CoinAbcCommonFactorization.hpp \
     553        CoinAbcCommon.hpp AbcDualRowDantzig.hpp AbcDualRowPivot.hpp \
     554        AbcDualRowSteepest.hpp CoinAbcFactorization.hpp \
     555        CoinAbcHelperFunctions.hpp AbcMatrix.hpp AbcCommon.hpp \
     556        AbcSimplexDual.hpp AbcSimplexPrimal.hpp AbcSimplex.hpp \
     557        CoinAbcBaseFactorization1.cpp CoinAbcBaseFactorization2.cpp \
     558        CoinAbcBaseFactorization3.cpp CoinAbcBaseFactorization4.cpp \
     559        CoinAbcBaseFactorization5.cpp CoinAbcDenseFactorization.cpp \
     560        CoinAbcDenseFactorization.hpp AbcDualRowDantzig.cpp \
     561        AbcDualRowPivot.cpp AbcDualRowSteepest.cpp \
     562        CoinAbcFactorization1.cpp CoinAbcFactorization2.cpp \
     563        CoinAbcFactorization3.cpp CoinAbcFactorization4.cpp \
     564        CoinAbcFactorization5.cpp CoinAbcHelperFunctions.cpp \
     565        AbcMatrix.cpp AbcSimplex.cpp AbcSimplexDual.cpp \
     566        AbcSimplexPrimal.cpp AbcSimplexParallel.cpp \
     567        AbcPrimalColumnDantzig.cpp AbcPrimalColumnDantzig.hpp \
     568        AbcPrimalColumnPivot.cpp AbcPrimalColumnPivot.hpp \
     569        AbcPrimalColumnSteepest.cpp AbcPrimalColumnSteepest.hpp \
     570        CoinAbcSmallFactorization1.cpp CoinAbcSmallFactorization2.cpp \
     571        CoinAbcSmallFactorization3.cpp CoinAbcSmallFactorization4.cpp \
     572        CoinAbcSmallFactorization5.cpp \
     573        CoinAbcOrderedFactorization1.cpp \
     574        CoinAbcOrderedFactorization2.cpp \
     575        CoinAbcOrderedFactorization3.cpp \
     576        CoinAbcOrderedFactorization4.cpp \
     577        CoinAbcOrderedFactorization5.cpp AbcSimplexFactorization.cpp \
     578        AbcSimplexFactorization.hpp AbcNonLinearCost.cpp \
     579        AbcNonLinearCost.hpp Idiot.cpp Idiot.hpp IdiSolve.cpp \
    504580        $(am__append_1) $(am__append_2) $(am__append_3) \
    505581        $(am__append_4) $(am__append_5)
     
    558634        ClpQuadraticObjective.hpp ClpSimplex.hpp \
    559635        ClpSimplexNonlinear.hpp ClpSimplexOther.hpp ClpSimplexDual.hpp \
    560         ClpSimplexPrimal.hpp ClpSolve.hpp CbcOrClpParam.hpp Idiot.hpp \
     636        ClpSimplexPrimal.hpp ClpSolve.hpp CbcOrClpParam.hpp \
     637        AbcSimplex.hpp CoinAbcCommon.hpp AbcCommon.hpp Idiot.hpp \
    561638        $(am__append_8) $(am__append_9) $(am__append_10) \
    562639        $(am__append_11) $(am__append_12) CbcOrClpParam.cpp
     
    689766        -rm -f *.tab.c
    690767
     768@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/AbcDualRowDantzig.Plo@am__quote@
     769@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/AbcDualRowPivot.Plo@am__quote@
     770@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/AbcDualRowSteepest.Plo@am__quote@
     771@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/AbcMatrix.Plo@am__quote@
     772@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/AbcNonLinearCost.Plo@am__quote@
     773@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/AbcPrimalColumnDantzig.Plo@am__quote@
     774@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/AbcPrimalColumnPivot.Plo@am__quote@
     775@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/AbcPrimalColumnSteepest.Plo@am__quote@
     776@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/AbcSimplex.Plo@am__quote@
     777@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/AbcSimplexDual.Plo@am__quote@
     778@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/AbcSimplexFactorization.Plo@am__quote@
     779@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/AbcSimplexParallel.Plo@am__quote@
     780@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/AbcSimplexPrimal.Plo@am__quote@
     781@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/AbcWarmStart.Plo@am__quote@
    691782@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/CbcOrClpParam.Po@am__quote@
    692783@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/ClpCholeskyBase.Plo@am__quote@
     
    739830@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/ClpSolve.Plo@am__quote@
    740831@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/Clp_C_Interface.Plo@am__quote@
     832@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/CoinAbcBaseFactorization1.Plo@am__quote@
     833@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/CoinAbcBaseFactorization2.Plo@am__quote@
     834@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/CoinAbcBaseFactorization3.Plo@am__quote@
     835@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/CoinAbcBaseFactorization4.Plo@am__quote@
     836@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/CoinAbcBaseFactorization5.Plo@am__quote@
     837@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/CoinAbcDenseFactorization.Plo@am__quote@
     838@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/CoinAbcFactorization1.Plo@am__quote@
     839@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/CoinAbcFactorization2.Plo@am__quote@
     840@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/CoinAbcFactorization3.Plo@am__quote@
     841@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/CoinAbcFactorization4.Plo@am__quote@
     842@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/CoinAbcFactorization5.Plo@am__quote@
     843@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/CoinAbcHelperFunctions.Plo@am__quote@
     844@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/CoinAbcOrderedFactorization1.Plo@am__quote@
     845@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/CoinAbcOrderedFactorization2.Plo@am__quote@
     846@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/CoinAbcOrderedFactorization3.Plo@am__quote@
     847@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/CoinAbcOrderedFactorization4.Plo@am__quote@
     848@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/CoinAbcOrderedFactorization5.Plo@am__quote@
     849@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/CoinAbcSmallFactorization1.Plo@am__quote@
     850@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/CoinAbcSmallFactorization2.Plo@am__quote@
     851@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/CoinAbcSmallFactorization3.Plo@am__quote@
     852@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/CoinAbcSmallFactorization4.Plo@am__quote@
     853@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/CoinAbcSmallFactorization5.Plo@am__quote@
    741854@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/IdiSolve.Plo@am__quote@
    742855@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/Idiot.Plo@am__quote@
  • trunk/Clp/src/unitTest.cpp

    r1726 r1878  
    2626#include "CoinFloatEqual.hpp"
    2727
     28#if CLP_HAS_ABC
     29#include "CoinAbcCommon.hpp"
     30#endif
     31#ifdef ABC_INHERIT
     32#include "CoinAbcFactorization.hpp"
     33#endif
    2834#include "ClpFactorization.hpp"
    2935#include "ClpSimplex.hpp"
     
    4450#include "ClpPresolve.hpp"
    4551#include "Idiot.hpp"
    46 
     52#if FACTORIZATION_STATISTICS
     53extern double ftranTwiddleFactor1X;
     54extern double ftranTwiddleFactor2X;
     55extern double ftranFTTwiddleFactor1X;
     56extern double ftranFTTwiddleFactor2X;
     57extern double btranTwiddleFactor1X;
     58extern double btranTwiddleFactor2X;
     59extern double ftranFullTwiddleFactor1X;
     60extern double ftranFullTwiddleFactor2X;
     61extern double btranFullTwiddleFactor1X;
     62extern double btranFullTwiddleFactor2X;
     63extern double denseThresholdX;
     64extern double twoThresholdX;
     65extern double minRowsSparse;
     66extern double largeRowsSparse;
     67extern double mediumRowsDivider;
     68extern double mediumRowsMinCount;
     69extern double largeRowsCount;
     70#endif
    4771
    4872//#############################################################################
     
    251275               << "        If specified, then netlib testset run as well as the nitTest.\n";
    252276}
    253 
     277#if FACTORIZATION_STATISTICS
     278int loSizeX=-1;
     279int hiSizeX=1000000;
     280#endif
    254281//----------------------------------------------------------------
     282#ifndef ABC_INHERIT
     283#define AnySimplex ClpSimplex
     284#else
     285#include "AbcSimplex.hpp"
     286#define AnySimplex AbcSimplex
     287#endif
    255288int mainTest (int argc, const char *argv[], int algorithm,
    256               ClpSimplex empty, ClpSolve solveOptionsIn,
     289              AnySimplex empty, ClpSolve solveOptionsIn,
    257290              int switchOffValue, bool doVector)
    258291{
     
    271304          }
    272305     }
    273 
     306     int numberFailures=0;
    274307     // define valid parameter keywords
    275308     std::set<std::string> definedKeyWords;
     
    318351     else
    319352          dirNetlib = dirsep == '/' ? "../../Data/Netlib/" : "..\\..\\Data\\Netlib\\";
     353#if 0 //FACTORIZATION_STATISTICS==0
    320354     if (!empty.numberRows()) {
    321355          testingMessage( "Testing ClpSimplex\n" );
    322356          ClpSimplexUnitTest(dirSample);
    323357     }
     358#endif
    324359     if (parms.find("-netlib") != parms.end() || empty.numberRows()) {
    325360          unsigned int m;
    326 
     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
    327391          // Define test problems:
    328392          //   mps names,
     
    352416          }
    353417          if (!empty.numberRows()) {
     418#if 1
    354419               mpsName.push_back("25fv47");
    355420               min.push_back(true);
     
    363428               nRows.push_back(2263);
    364429               nCols.push_back(9799);
    365                objValueTol.push_back(1.e-10);
     430               objValueTol.push_back(1.e-8);
    366431               objValue.push_back(9.8722419241E+05);
    367432               bestStrategy.push_back(3);
     
    370435               nRows.push_back(75);
    371436               nCols.push_back(83);
    372                objValueTol.push_back(1.e-10);
     437               objValueTol.push_back(1.e-8);
    373438               objValue.push_back(-3.0812149846e+01);
    374439               bestStrategy.push_back(3);
     
    401466               objValue.push_back(-2.5811392641e+03);
    402467               bestStrategy.push_back(3);
     468#endif
    403469               mpsName.push_back("pilot87");
    404470               min.push_back(true);
     
    408474               objValue.push_back(3.0171072827e+02);
    409475               bestStrategy.push_back(0);
     476#if 1
    410477               mpsName.push_back("adlittle");
    411478               min.push_back(true);
     
    9871054               objValue.push_back(1.3044763331E+00);
    9881055               bestStrategy.push_back(3);
     1056#endif
    9891057          } else {
    9901058               // Just testing one
     
    10151083          // Loop once for each Mps File
    10161084          for (m = 0; m < mpsName.size(); m++ ) {
     1085#if FACTORIZATION_STATISTICS
     1086               if (nRows[m]<loSizeX||nRows[m]>=hiSizeX) {
     1087                 std::cerr << "  skipping mps file: " << mpsName[m]
     1088                           <<" as "<<nRows[m]
     1089                         << " (" << m + 1 << " out of " << mpsName.size() << ")" << std::endl;
     1090                 continue;
     1091               }
     1092#endif
    10171093               std::cerr << "  processing mps file: " << mpsName[m]
    10181094                         << " (" << m + 1 << " out of " << mpsName.size() << ")" << std::endl;
    1019 
    1020                ClpSimplex solutionBase = empty;
     1095               AnySimplex solutionBase = empty;
    10211096               std::string fn = dirNetlib + mpsName[m];
    10221097               if (!empty.numberRows() || algorithm < 6) {
     
    10601135                              double time2 = CoinCpuTime() - time1;
    10611136                              testTime[iTest] = time2;
    1062                               printf("Took %g seconds - status %d\n", time2, solution.problemStatus());
     1137                              printf("Finished %s Took %g seconds (%d iterations) - status %d\n",
     1138                                     mpsName[m].c_str(),time2, solution.problemStatus(),solution.numberIterations());
    10631139                              if (solution.problemStatus())
    10641140                                   testTime[iTest] = 1.0e20;
     
    10891165               }
    10901166               double time1 = CoinCpuTime();
    1091                ClpSimplex solution = solutionBase;
     1167               AnySimplex solution = solutionBase;
    10921168#if 0
    10931169               solution.setOptimizationDirection(-1);
     
    11511227                    }
    11521228               }
     1229#if FACTORIZATION_STATISTICS
     1230               double timesOne[3*3*3*3*3];
     1231               memset(timesOne,0,sizeof(timesOne));
     1232               int iterationsOne[3*3*3*3*3];
     1233               memset(iterationsOne,0,sizeof(iterationsOne));
     1234               AnySimplex saveModel(solution);
     1235               double time2;
     1236#if 0
     1237               solution.initialSolve(solveOptions);
     1238               time2 = CoinCpuTime() - time1;
     1239               timeTaken += time2;
     1240               printf("%s took %g seconds using algorithm %s\n", fn.c_str(), time2, nameAlgorithm.c_str());
     1241#endif
     1242#define loA 1
     1243#define loB 0
     1244#define loC 0
     1245#define loD 1
     1246#define loE 1
     1247#define hiA 2
     1248#define hiB 1
     1249#define hiC 1
     1250#define hiD 2
     1251#define hiE 2
     1252               time2 = CoinCpuTime();
     1253               for (int iA=loA;iA<hiA;iA++) {
     1254                 *externalA=useA[iA];
     1255                 for (int iB=loB;iB<hiB;iB++) {
     1256                   *externalB=useB[iB];
     1257                   for (int iC=loC;iC<hiC;iC++) {
     1258                     *externalC=useC[iC];
     1259                     for (int iD=loD;iD<hiD;iD++) {
     1260                       *externalD=useD[iD];
     1261                       for (int iE=loE;iE<hiE;iE++) {
     1262                         *externalE=useE[iE];
     1263                         solution=saveModel;
     1264                         solution.initialSolve(solveOptions);
     1265                         double time3=CoinCpuTime();
     1266                         printf("Finished %s Took %g seconds (%d iterations) - status %d\n",
     1267                                mpsName[m].c_str(),time3-time2,solution.numberIterations(), solution.problemStatus());
     1268                         iterationsOne[iA+3*iB+9*iC+27*iD+81*iE]=solution.numberIterations();
     1269                         timesOne[iA+3*iB+9*iC+27*iD+81*iE]=time3-time2;
     1270                         times[iA+3*iB+9*iC+27*iD+81*iE]+=time3-time2;
     1271                         time2=time3;
     1272                       }
     1273                     }
     1274                   }
     1275                 }
     1276               }
     1277               double bestTime=1.0e100;
     1278               int iBestTime=-1;
     1279               double bestTimePer=1.0e100;
     1280               int iBestTimePer=-1;
     1281               int bestIterations=1000000;
     1282               int iBestIterations=-1;
     1283               printf("TTimes ");
     1284               for (int i=0;i<3*3*3*3*3;i++) {
     1285                 // skip if not done
     1286                 if (!iterationsOne[i])
     1287                   continue;
     1288                 double average=timesOne[i]/static_cast<double>(iterationsOne[i]);
     1289                 if (timesOne[i]<bestTime) {
     1290                   bestTime=timesOne[i];
     1291                   iBestTime=i;
     1292                 }
     1293                 if (average<bestTimePer) {
     1294                   bestTimePer=average;
     1295                   iBestTimePer=i;
     1296                 }
     1297                 if (iterationsOne[i]<bestIterations) {
     1298                   bestIterations=iterationsOne[i];
     1299                   iBestIterations=i;
     1300                 }
     1301                 printf("%.2f ",timesOne[i]);
     1302               }
     1303               printf("\n");
     1304               int iA,iB,iC,iD,iE;
     1305               iA=iBestTime;
     1306               iE=iA/81;
     1307               iA-=81*iE;
     1308               iD=iA/27;
     1309               iA-=27*iD;
     1310               iC=iA/9;
     1311               iA-=9*iC;
     1312               iB=iA/3;
     1313               iA-=3*iB;
     1314               printf("Best time %.2f %s=%g %s=%g %s=%g %s=%g %s=%g\n",bestTime,choiceA,useA[iA],
     1315                      choiceB,useB[iB],choiceC,useC[iC],choiceD,useD[iD],choiceE,useE[iE]);
     1316               iA=iBestTimePer;
     1317               iE=iA/81;
     1318               iA-=81*iE;
     1319               iD=iA/27;
     1320               iA-=27*iD;
     1321               iC=iA/9;
     1322               iA-=9*iC;
     1323               iB=iA/3;
     1324               iA-=3*iB;
     1325               printf("Best time per iteration %g %s=%g %s=%g %s=%g %s=%g %s=%g\n",bestTimePer,choiceA,useA[iA],
     1326                      choiceB,useB[iB],choiceC,useC[iC],choiceD,useD[iD],choiceE,useE[iE]);
     1327               iA=iBestIterations;
     1328               iE=iA/81;
     1329               iA-=81*iE;
     1330               iD=iA/27;
     1331               iA-=27*iD;
     1332               iC=iA/9;
     1333               iA-=9*iC;
     1334               iB=iA/3;
     1335               iA-=3*iB;
     1336               printf("Best iterations %d %s=%g %s=%g %s=%g %s=%g %s=%g\n",bestIterations,choiceA,useA[iA],
     1337                      choiceB,useB[iB],choiceC,useC[iC],choiceD,useD[iD],choiceE,useE[iE]);
     1338#else
    11531339               solution.initialSolve(solveOptions);
    11541340               double time2 = CoinCpuTime() - time1;
    11551341               timeTaken += time2;
    11561342               printf("%s took %g seconds using algorithm %s\n", fn.c_str(), time2, nameAlgorithm.c_str());
     1343#endif
    11571344               // Test objective solution value
    11581345               {
     1346                    if (!solution.isProvenOptimal()) {
     1347                      std::cerr <<"** NOT OPTIMAL ";
     1348                      numberFailures++;
     1349                    }
    11591350                    double soln = solution.objectiveValue();
    11601351                    CoinRelFltEq eq(objValueTol[m]);
    11611352                    std::cerr << soln << ",  " << objValue[m] << " diff " <<
    11621353                              soln - objValue[m] << std::endl;
    1163                     if(!eq(soln, objValue[m]))
     1354                    if(!eq(soln, objValue[m])) {
    11641355                         printf("** difference fails\n");
     1356                         numberFailures++;
     1357                    }
    11651358               }
    11661359          }
    11671360          printf("Total time %g seconds\n", timeTaken);
     1361#if FACTORIZATION_STATISTICS
     1362          double bestTime=1.0e100;
     1363          int iBestTime=-1;
     1364          for (int i=0;i<3*3*3*3*3;i++) {
     1365            if (times[i]&&times[i]<bestTime) {
     1366              bestTime=times[i];
     1367              iBestTime=i;
     1368            }
     1369          }
     1370          for (int iA=0;iA<hiA;iA++) {
     1371            for (int iB=0;iB<hiB;iB++) {
     1372              for (int iC=0;iC<hiC;iC++) {
     1373                for (int iD=loD;iD<hiD;iD++) {
     1374                  for (int iE=loE;iE<hiE;iE++) {
     1375                    int k=iA+3*iB+9*iC+27*iD+81*iE;
     1376                    double thisTime=times[k];
     1377                    if (thisTime) {
     1378                      if (k==iBestTime)
     1379                        printf("** ");
     1380                      printf("%s=%g %s=%g %s=%g %s=%g %s=%g - time %.2f\n",
     1381                             choiceA,useA[iA],
     1382                             choiceB,useB[iB],choiceC,useC[iC],
     1383                             choiceD,useD[iD],choiceE,useE[iE],thisTime);
     1384                    }
     1385                  }
     1386                }
     1387              }
     1388            }
     1389          }
     1390          //exit(0);
     1391#endif
    11681392     } else {
    11691393          testingMessage( "***Skipped Testing on netlib    ***\n" );
    11701394          testingMessage( "***use -netlib to test class***\n" );
    11711395     }
    1172 
    1173      testingMessage( "All tests completed successfully\n" );
     1396     if (!numberFailures)
     1397       testingMessage( "All tests completed successfully\n" );
     1398     else
     1399       testingMessage( "Some tests failed\n" );
    11741400     return 0;
    11751401}
Note: See TracChangeset for help on using the changeset viewer.