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

minor changes to implement Aboca

File:
1 edited

Legend:

Unmodified
Added
Removed
  • 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;
Note: See TracChangeset for help on using the changeset viewer.