Changeset 2070 for trunk


Ignore:
Timestamp:
Nov 18, 2014 6:12:54 AM (4 years ago)
Author:
forrest
Message:

for some cbc ideas

Location:
trunk/Clp/src
Files:
20 edited

Legend:

Unmodified
Added
Removed
  • trunk/Clp/src/AbcSimplex.hpp

    r2026 r2070  
    2121#include "ClpSolve.hpp"
    2222#include "CoinAbcCommon.hpp"
    23 class ClpSimplex;
     23#include "ClpSimplex.hpp"
    2424class AbcDualRowPivot;
    2525class AbcPrimalColumnPivot;
     
    6060#include <pthread.h>
    6161#endif
    62 typedef struct {
    63   double result;
    64   //const CoinIndexedVector * constVector; // can get rid of
    65   //CoinIndexedVector * vectors[2]; // can get rid of
    66   void * extraInfo;
    67   int status;
    68   int stuff[4];
    69 } CoinAbcThreadInfo;
    70 #include "ClpSimplex.hpp"
    7162class AbcSimplex : public ClpSimplex {
    7263  friend void AbcSimplexUnitTest(const std::string & mpsDir);
     
    869860  /// Move status and solution across
    870861  void moveInfo(const AbcSimplex & rhs, bool justStatus = false);
     862#ifndef NUMBER_THREADS
    871863#define NUMBER_THREADS 3
     864#endif
    872865#if ABC_PARALLEL==1
    873866  // For waking up thread
     
    878871  inline int whichLocked(int thread=0) const
    879872  { return locked_[thread];}
    880   inline CoinAbcThreadInfo * threadInfoPointer(int thread=0)
     873  inline CoinThreadInfo * threadInfoPointer(int thread=0)
    881874  { return threadInfo_+thread;}
    882875  void startParallelStuff(int type);
     
    885878  int whichThread() const;
    886879#elif ABC_PARALLEL==2
    887   //inline CoinAbcThreadInfo * threadInfoPointer(int thread=0)
     880  //inline CoinThreadInfo * threadInfoPointer(int thread=0)
    888881  //{ return threadInfo_+thread;}
    889882#endif
     
    12701263  pthread_mutex_t mutex_[3*NUMBER_THREADS];
    12711264  pthread_barrier_t barrier_;
    1272   CoinAbcThreadInfo threadInfo_[NUMBER_THREADS];
     1265  CoinThreadInfo threadInfo_[NUMBER_THREADS];
    12731266  pthread_t abcThread_[NUMBER_THREADS];
    12741267  int locked_[NUMBER_THREADS];
    12751268  int stopStart_;
    12761269#elif ABC_PARALLEL==2
    1277   //CoinAbcThreadInfo threadInfo_[NUMBER_THREADS];
     1270  //CoinThreadInfo threadInfo_[NUMBER_THREADS];
    12781271#endif
    12791272  //@}
  • trunk/Clp/src/AbcSimplexDual.cpp

    r2042 r2070  
    23862386#else
    23872387        whichArray[0]=1;
    2388         CoinAbcThreadInfo info;
     2388        CoinThreadInfo info;
    23892389        info.status=1;
    23902390        info.stuff[0]=whichArray[1];
     
    57705770  AbcSimplexDual * dual = reinterpret_cast<AbcSimplexDual *>(simplex);
    57715771  int whichThread=dual->whichThread();
    5772   CoinAbcThreadInfo * threadInfo = dual->threadInfoPointer(whichThread);
     5772  CoinThreadInfo * threadInfo = dual->threadInfoPointer(whichThread);
    57735773  pthread_mutex_lock(dual->mutexPointer(2,whichThread));
    57745774  pthread_barrier_wait(dual->barrierPointer());
  • trunk/Clp/src/AbcSimplexParallel.cpp

    r1910 r2070  
    16361636#if ABC_PARALLEL==1
    16371637  // redo so can pass info in with void *
    1638   CoinAbcThreadInfo * infoP = dual->threadInfoPointer();
     1638  CoinThreadInfo * infoP = dual->threadInfoPointer();
    16391639  int cpuMask=((1<<dual->parallelMode())-1);
    16401640  cpuMask += cpuMask<<5;
    16411641  dual->setStopStart(cpuMask);
    16421642#endif
    1643   CoinAbcThreadInfo info[NUMBER_BLOCKS];
     1643  CoinThreadInfo info[NUMBER_BLOCKS];
    16441644  const int * starts;
    16451645  if (useRowCopy)
  • trunk/Clp/src/CbcOrClpParam.cpp

    r2069 r2070  
    964964          returnCode = 1;
    965965     } else {
     966          printArray[0] = '\0';
     967          if (value==intValue_)
     968            return printArray;
    966969          int oldValue = intValue_;
    967970          intValue_ = value;
     
    11251128     return printArray;
    11261129}
     1130// Sets current parameter option using string with message
     1131const char *
     1132CbcOrClpParam::setCurrentOptionWithMessage ( const std::string value )
     1133{
     1134     int action = parameterOption(value);
     1135     char current[100];
     1136     printArray[0] = '\0';
     1137     if (action >= 0) {
     1138         if (action == currentKeyWord_)
     1139           return NULL;
     1140         if (currentKeyWord_>=0&&(fakeKeyWord_<=0||currentKeyWord_<fakeKeyWord_))
     1141           strcpy(current,definedKeyWords_[currentKeyWord_].c_str());
     1142         else if (currentKeyWord_<0)
     1143           sprintf(current,"minus%d",-currentKeyWord_-1000);
     1144         else
     1145           sprintf(current,"plus%d",currentKeyWord_-1000);
     1146         sprintf(printArray, "Option for %s changed from %s to %s",
     1147                 name_.c_str(), current, value.c_str());
     1148          currentKeyWord_ = action;
     1149     } else {
     1150         sprintf(printArray, "Option for %s given illegal value %s",
     1151                 name_.c_str(), value.c_str());
     1152     }
     1153     return printArray;
     1154}
    11271155void
    11281156CbcOrClpParam::setIntValue ( int value )
     
    11361164     }
    11371165}
     1166const char *
     1167CbcOrClpParam::setIntValueWithMessage ( int value )
     1168{
     1169     printArray[0] = '\0';
     1170     if (value < lowerIntValue_ || value > upperIntValue_) {
     1171          sprintf(printArray, "%d was provided for %s - valid range is %d to %d",
     1172               value,name_.c_str(),lowerIntValue_,upperIntValue_);
     1173     } else {
     1174         if (value==intValue_)
     1175           return NULL;
     1176          sprintf(printArray, "%s was changed from %d to %d",
     1177                  name_.c_str(), intValue_, value);
     1178          intValue_ = value;
     1179     }
     1180     return printArray;
     1181}
    11381182void
    11391183CbcOrClpParam::setDoubleValue ( double value )
     
    11461190          doubleValue_ = value;
    11471191     }
     1192}
     1193const char *
     1194CbcOrClpParam::setDoubleValueWithMessage ( double value )
     1195{
     1196     printArray[0] = '\0';
     1197     if (value < lowerDoubleValue_ || value > upperDoubleValue_) {
     1198          sprintf(printArray, "%g was provided for %s - valid range is %g to %g",
     1199               value,name_.c_str(),lowerDoubleValue_,upperDoubleValue_);
     1200     } else {
     1201         if (value==doubleValue_)
     1202           return NULL;
     1203          sprintf(printArray, "%s was changed from %g to %g",
     1204                  name_.c_str(), doubleValue_, value);
     1205          doubleValue_ = value;
     1206     }
     1207     return printArray;
    11481208}
    11491209void
     
    16691729     parameters[numberParameters-1].append("variable");
    16701730     parameters[numberParameters-1].append("forcevariable");
     1731     parameters[numberParameters-1].append("conflict");
    16711732     parameters[numberParameters-1].setLonghelp
    16721733     (
     
    19442005         \n\t7 run up to 2 times if solution found 4 otherwise, \
    19452006         \n\t>10 All only at root (DivingC normal as value-10), \
    1946          \n\t>10 All with value-20)."
     2007         \n\t>20 All with value-20)."
    19472008     );
    19482009     parameters[numberParameters-1].setIntValue(-1);
     
    19562017         \n\t1-3 allow fixing of satisfied integers (but not at bound) \
    19572018         \n\t1 switch off above for that dive if goes infeasible \
    1958          \n\t2 switch off above permanenty if goes infeasible"
     2019         \n\t2 switch off above permanently if goes infeasible"
    19592020     );
    19602021     parameters[numberParameters-1].setIntValue(100);
     
    24942555     parameters[numberParameters-1].append("endonlyroot");
    24952556     parameters[numberParameters-1].append("endcleanroot");
    2496      parameters[numberParameters-1].append("endbothroot");
     2557     parameters[numberParameters-1].append("root");
    24972558     parameters[numberParameters-1].append("endonly");
    24982559     parameters[numberParameters-1].append("endclean");
  • trunk/Clp/src/CbcOrClpParam.hpp

    r2069 r2070  
    386386     /// Sets current parameter option using string
    387387     void setCurrentOption (const std::string value );
     388     /// Sets current parameter option using string with message
     389     const char * setCurrentOptionWithMessage (const std::string value );
    388390     /// Returns current parameter option position
    389391     int currentOptionAsInteger (  ) const ;
     
    395397     /// Sets int value
    396398     void setIntValue ( int value );
     399     /// Sets int value with message
     400     const char * setIntValueWithMessage ( int value );
    397401     inline int intValue () const {
    398402          return intValue_;
     
    400404     /// Sets double value
    401405     void setDoubleValue ( double value );
     406     /// Sets double value with message
     407     const char * setDoubleValueWithMessage ( double value );
    402408     inline double doubleValue () const {
    403409          return doubleValue_;
  • trunk/Clp/src/ClpDualRowPivot.hpp

    r1732 r2070  
    5555            (e.g. check for infeasible)
    5656         4) as 2 but restore weights from previous snapshot
    57          5) for strong branching - initialize  , infeasibilities
     57         5) for strong branching - initialize to 1 , infeasibilities
     58         6) scale back
     59         7) for strong branching - initialize full weights , infeasibilities
    5860     */
    5961     virtual void saveWeights(ClpSimplex * model, int mode);
  • trunk/Clp/src/ClpDualRowSteepest.cpp

    r1732 r2070  
    804804   2) after factorization
    805805   3) just redo infeasibilities
    806    4) restore weights
     806   4) as 2 but restore weights from previous snapshot
     807   5) for strong branching - initialize to 1 , infeasibilities
     808   6) scale back
     809   7) for strong branching - initialize full weights , infeasibilities
    807810*/
    808811void
     
    827830                    }
    828831                    state_ = 1;
     832                    // clear marker on savedWeights_
     833                    assert (savedWeights_);
     834                    if (savedWeights_->packedMode()) {
     835                      savedWeights_->setPackedMode(false);
     836                      savedWeights_->setNumElements(0);
     837                    }
    829838               } else {
    830839                    // size has changed - clear everything
     
    844853     } else if (mode == 2 || mode == 4 || mode >= 5) {
    845854          // restore
    846           if (!weights_ || state_ == -1 || mode == 5) {
     855          if (!weights_ || state_ == -1 || mode == 5 || mode == 7) {
    847856               // initialize weights
    848857               delete [] weights_;
    849858               delete alternateWeights_;
    850859               weights_ = new double[numberRows];
     860               // initialize to 1.0 (can we do better?)
     861               for (i = 0; i < numberRows; i++) {
     862                 weights_[i] = 1.0;
     863               }
    851864               alternateWeights_ = new CoinIndexedVector();
    852865               // enough space so can use it for factorization
     
    854867                                          model_->factorization()->maximumPivots());
    855868               if (mode_ != 1 || mode == 5) {
    856                     // initialize to 1.0 (can we do better?)
    857                     for (i = 0; i < numberRows; i++) {
    858                          weights_[i] = 1.0;
    859                     }
    860869               } else {
    861870                    CoinIndexedVector * temp = new CoinIndexedVector();
     
    864873                    double * array = alternateWeights_->denseVector();
    865874                    int * which = alternateWeights_->getIndices();
    866                     for (i = 0; i < numberRows; i++) {
     875                    int firstRow=0;
     876                    int lastRow=numberRows;
     877                    if (mode==7) {
     878                      // use info passed in
     879                      firstRow=model->spareIntArray_[0];
     880                      lastRow=model->spareIntArray_[1];
     881                    }
     882                    for (i = firstRow; i < lastRow; i++) {
    867883                         double value = 0.0;
    868884                         array[0] = 1.0;
     
    886902               savedWeights_ = new CoinIndexedVector();
    887903               savedWeights_->reserve(numberRows);
     904               for (int i=0;i<model_->numberRows();i++)
     905                 savedWeights_->denseVector()[i]=1.0;
    888906
    889907               double * array = savedWeights_->denseVector();
    890908               int * which = savedWeights_->getIndices();
    891                for (i = 0; i < numberRows; i++) {
     909               for (int i = 0; i < numberRows; i++) {
    892910                    array[i] = weights_[i];
    893911                    which[i] = pivotVariable[i];
    894912               }
     913               if (mode==7) {
     914                 savedWeights_->setNumElements(numberRows); //flag as special
     915                 savedWeights_->setPackedMode(true);
     916               }
    895917          } else if (mode != 6) {
    896918               int * which = alternateWeights_->getIndices();
     
    919941                    back[iSeq] = i;
    920942               }
    921                for (i = 0; i < numberRows; i++) {
     943               int firstRow=0;
     944               int lastRow=numberRows;
     945               if (mode==7) {
     946                 // use info passed in
     947                 firstRow=model->spareIntArray_[0];
     948                 lastRow=model->spareIntArray_[1];
     949               }
     950               for (i = firstRow; i < lastRow; i++) {
    922951                    int iPivot = pivotVariable[i];
    923952                    iPivot = back[iPivot];
     
    929958                         // odd
    930959                         weights_[i] = 1.0;
     960                         //printf("bad pivot row %d (%d->%d) iPivot %d\n",i,firstRow,lastRow,iPivot);
    931961                    }
    932962               }
     963#if 0
     964               printf("mode %d mode_ %d state_ %d\n",mode,mode_,state_);
     965               if (!model_->numberIterations()) {
     966                 for (int k=0;k<numberRows;k+=10)
     967                   printf("%d %g ",k,weights_[k]);
     968                 printf("\n");
     969               }
     970#endif
    933971          } else {
    934972               // mode 6 - scale back weights as primal errors
     
    10031041          }
    10041042     }
     1043     // see where coming from
     1044     if (mode==2&&!model->numberIterations()) {
     1045       int options=model->specialOptions();
     1046       if ((options&16384)!=0&&true) {
     1047         // fast of some sort - could be clever???
     1048         // for now initialize
     1049         if ((options&524288)!=0&&false) {
     1050           // fathom
     1051           for (int i = 0; i < numberRows; i++)
     1052             weights_[i] = 1.0;
     1053         } else if (true) {
     1054           // strong
     1055           for (int i = 0; i < numberRows; i++)
     1056             weights_[i] = 1.0;
     1057         }
     1058       }
     1059     }
     1060}
     1061// Pass in saved weights
     1062void
     1063ClpDualRowSteepest::passInSavedWeights(const CoinIndexedVector * saved)
     1064{
     1065  delete savedWeights_;
     1066  savedWeights_=new CoinIndexedVector(*saved);
    10051067}
    10061068// Gets rid of last update
  • trunk/Clp/src/ClpDualRowSteepest.hpp

    r1665 r2070  
    5555     */
    5656     virtual void saveWeights(ClpSimplex * model, int mode);
     57     /// Pass in saved weights
     58     void passInSavedWeights(const CoinIndexedVector * saved);
     59     /// Get saved weights
     60     inline CoinIndexedVector * savedWeights()
     61     { return savedWeights_;}
    5762     /// Gets rid of last update
    5863     virtual void unrollWeights();
     
    105110          return mode_;
    106111     }
     112     /// Set mode
     113     inline void setMode(int mode) {
     114          mode_ = mode;
     115     }
    107116     /// Set/ get persistence
    108117     inline void setPersistence(Persistence life) {
  • trunk/Clp/src/ClpPresolve.cpp

    r1977 r2070  
    10841084                              break;
    10851085                    }
     1086                    if (zerocost) {
     1087                      possibleBreak;
     1088                         paction_ = do_tighten_action::presolve(prob, paction_);
     1089                         if (prob->status_)
     1090                              break;
     1091                         printProgress('J',iLoop+1);
     1092                    }
    10861093                    if (dual && whichPass == 1) {
    10871094                         // this can also make E rows so do one bit here
     
    11081115                    }
    11091116
    1110                     if (zerocost) {
    1111                       possibleBreak;
    1112                          paction_ = do_tighten_action::presolve(prob, paction_);
    1113                          if (prob->status_)
    1114                               break;
    1115                          printProgress('J',iLoop+1);
    1116                     }
    11171117#ifndef NO_FORCING
    11181118                    if (forcing) {
  • trunk/Clp/src/ClpSimplex.cpp

    r2030 r2070  
    43604360     }
    43614361#endif
     4362     if ((moreSpecialOptions_&4194304)!=0) {
     4363       // preset tolerances were changed
     4364       moreSpecialOptions_ &= ~4194304;
     4365       primalTolerance_=1.0e-7;
     4366       dblParam_[ClpPrimalTolerance]=1.0e-7;
     4367       dualTolerance_=1.0e-7;
     4368       dblParam_[ClpDualTolerance]=1.0e-7;
     4369     }
    43624370     // ray may be null if in branch and bound
    43634371     if (rowScale_) {
     
    54705478          moreSpecialOptions_ &= ~256;
    54715479          baseIteration_ = 0;
     5480          bool inCbcOrOther = (specialOptions_ & 0x03000000) != 0;
     5481          if (inCbcOrOther) {
     5482            delete [] ray_;
     5483            ray_ = NULL;
     5484          }
    54725485          if (saveObjective != objective_) {
    54735486               // We changed objective to see if infeasible
     
    55595572     //if (problemStatus_==1&&lastAlgorithm==1)
    55605573     //returnCode=10; // so will do primal after postsolve
    5561      if (!problemStatus_) {
    5562           //assert (!numberPrimalInfeasibilities_);
    5563           //if (returnCode!=10)
    5564           //assert (!numberDualInfeasibilities_);
    5565      }
    55665574     return returnCode;
    55675575}
     
    1030110309          const int freq0 = 50;
    1030210310          int frequency;
    10303 #ifndef ABC_INHERIT
     10311#if ABC_CLP_DEFAULTS
    1030410312          const int cutoff2 = 100000;
    1030510313          const int freq1 = 200;
     
    1049710505     moreSpecialOptions_ |= 8;
    1049810506     int saveMaxIterations = maximumIterations();
    10499      setMaximumIterations((((moreSpecialOptions_&2048)==0) ? 100 : 2000)
    10500                           + 5 * (numberRows_ + numberColumns_));
     10507     setMaximumIterations((((moreSpecialOptions_&2048)==0) ? 200 : 2000)
     10508                          + 2 * (numberRows_ + numberColumns_));
    1050110509     double saveObjLimit;
    1050210510     getDblParam(ClpDualObjectiveLimit, saveObjLimit);
     
    1081310821     bool backtrack = false;
    1081410822     bool printing = handler_->logLevel() > 0;
     10823     // Say can stop without cleaning up in primal
     10824     moreSpecialOptions_ |= 2097152;
    1081510825     while (depth >= 0) {
    1081610826          // If backtrack get to correct depth
     
    1088610896          // Get fake bounds correctly
    1088710897          (static_cast<ClpSimplexDual *>(this))->changeBounds(3, NULL, dummyChange);
    10888           fastDual2(info);
     10898          int statusWeights = fastDual2(info);
     10899#if 0
     10900          if (statusWeights==100)
     10901            printf("weights ok\n");
     10902          else
     10903            printf("weights not ok\n");
     10904#endif
    1088910905#if 0
    1089010906          {
     
    1090510921#endif
    1090610922          // give up if odd
     10923          numberNodes++;
     10924          numberIterations += numberIterations_;
    1090710925          if (problemStatus_ > 1) {
    1090810926               info->nNodes_ = -1;
     
    1091310931               //abort();
    1091410932#endif
     10933               printf("too long in one solve (%d iterations) average %d iterations %d nodes\n",
     10934                      numberIterations_,numberIterations,numberNodes);
    1091510935               break;
    1091610936          }
    10917           numberNodes++;
    10918           numberIterations += numberIterations_;
     10937          if (numberNodes>20) {
     10938            if (numberIterations>200*numberNodes) {
     10939              info->nNodes_ = -2;
     10940              printf("too long average %d iterations %d nodes\n",
     10941                     numberIterations,numberNodes);
     10942              break;
     10943            }
     10944          } else {
     10945            if (numberIterations>4000) {
     10946              info->nNodes_ = -2;
     10947              printf("too long average %d iterations %d nodes\n",
     10948                     numberIterations,numberNodes);
     10949              break;
     10950            }
     10951          }
    1091910952          if ((numberNodes % 1000) == 0) {
    1092010953#ifdef COIN_DEVELOP
     
    1093610969#endif
    1093710970               if ((numberIterations*(numberRows_ + numberColumns_) > 5.0e10 ||
    10938                     numberNodes > 2.0e4) &&
     10971                    numberNodes > 1.5e4) &&
    1093910972                   (moreSpecialOptions_&4096)==0) {
    1094010973                    // give up
    10941                     info->nNodes_ = -1;
     10974                    info->nNodes_ = -10;
    1094210975#ifdef COIN_DEVELOP
    1094310976                    printf("OUCH giving up on nodes %d %d\n", numberNodes, numberIterations);
     
    1139211425     bool backtrack = false;
    1139311426     bool printing = handler_->logLevel() > 0;
     11427     // Say can stop without cleaning up in primal
     11428     moreSpecialOptions_ |= 2097152;
    1139411429#ifdef CHECK_PATH
    1139511430     if (startOptimal)
     
    1189011925     //int saveNumberFake = numberFake_;
    1189111926     int status = static_cast<ClpSimplexDual *> (this)->fastDual(true);
     11927     bool goodWeights=true;
    1189211928     //numberFake_ = saveNumberFake;
    1189311929     specialOptions_ &= ~524288; // say dont use solution
     
    1190111937               problemStatus_ = 0;
    1190211938          }
     11939     } else if (problemStatus_==10 &&
     11940                (moreSpecialOptions_&2097152) != 0) {
     11941          // not finished - may want to use anyway
     11942          checkPrimalSolution(rowActivityWork_, columnActivityWork_);
     11943          double limit = 0.0;
     11944          getDblParam(ClpDualObjectiveLimit, limit);
     11945          if (!numberPrimalInfeasibilities_ && objectiveValue()*optimizationDirection_ < limit) {
     11946               problemStatus_ = 11;
     11947          }
    1190311948     }
    1190411949     if (problemStatus_ == 10) {
     
    1190711952          //printf("Cleaning up with primal\n");
    1190811953          //lastAlgorithm=1;
     11954          goodWeights=false;
    1190911955          int savePerturbation = perturbation_;
    1191011956          int saveLog = handler_->logLevel();
     
    1195912005               perturbation_ = savePerturbation;
    1196012006               baseIteration_ = numberIterations_;
     12007               goodWeights=false;
    1196112008               status = static_cast<ClpSimplexPrimal *> (this)->primal(0);
    1196212009               baseIteration_ = 0;
     
    1198912036     }
    1199012037     status = problemStatus_;
    11991      if (!problemStatus_) {
     12038     if (!problemStatus_||problemStatus_==11) {
    1199212039          int j;
    1199312040          // Move solution to external array
     
    1200012047          if ((info->solverOptions_ & 1) != 0) {
    1200112048               // reduced costs
    12002                if (!columnScale_) {
    12003                     CoinMemcpyN(dj_, numberColumns_, reducedCost_);
    12004                } else {
    12005                     for (j = 0; j < numberColumns_; j++)
    12006                          reducedCost_[j] = dj_[j] * columnScale_[j+numberColumns_];
    12007                }
     12049               if (!problemStatus_) {
     12050                 if (!columnScale_) {
     12051                   CoinMemcpyN(dj_, numberColumns_, reducedCost_);
     12052                 } else {
     12053                   for (j = 0; j < numberColumns_; j++)
     12054                     reducedCost_[j] = dj_[j] * columnScale_[j+numberColumns_];
     12055                 }
     12056               } else {
     12057                 // zero djs
     12058                 memset(reducedCost_,0,numberColumns_*sizeof(double));
     12059                 problemStatus_=0;
     12060               }
    1200812061          }
    1200912062          if ((info->solverOptions_ & 2) != 0) {
     
    1203412087     CoinMemcpyN(save, numberTotal, upper_);
    1203512088#endif
     12089     if (goodWeights)
     12090       status=100;
    1203612091     return status;
    1203712092}
  • trunk/Clp/src/ClpSimplex.hpp

    r2025 r2070  
    2929class ClpDisasterHandler;
    3030class ClpConstraint;
     31/*
     32  May want to use Clp defaults so that with ABC defined but not used
     33  it behaves as Clp (and ABC used will be different than if not defined)
     34 */
     35#ifdef ABC_INHERIT
     36#ifndef ABC_CLP_DEFAULTS
     37#define ABC_CLP_DEFAULTS 0
     38#endif
     39#else
     40#undef ABC_CLP_DEFAULTS
     41#define ABC_CLP_DEFAULTS 1
     42#endif
    3143#ifdef CLP_HAS_ABC
    3244#include "AbcCommon.hpp"
     
    12221234         524288 bit - stop when primal feasible
    12231235         1048576 bit - don't perturb even if long time
     1236         2097152 bit - no primal in fastDual2 if feasible
     1237         4194304 bit - tolerances have been changed by code
    12241238     */
    12251239     inline void setMoreSpecialOptions(int value) {
     
    17101724#define DEVEX_TRY_NORM 1.0e-4
    17111725#define DEVEX_ADD_ONE 1.0
     1726#if defined(ABC_INHERIT) || defined(CBC_THREAD)
     1727// Use pthreads
     1728#include <pthread.h>
     1729typedef struct {
     1730  double result;
     1731  //const CoinIndexedVector * constVector; // can get rid of
     1732  //CoinIndexedVector * vectors[2]; // can get rid of
     1733  void * extraInfo;
     1734  void * extraInfo2;
     1735  int status;
     1736  int stuff[4];
     1737} CoinThreadInfo;
     1738class CoinPthreadStuff {
     1739public:
     1740  /**@name Constructors and destructor and copy */
     1741  //@{
     1742  /** Main constructor
     1743  */
     1744  CoinPthreadStuff (int numberThreads=0,
     1745                    void * parallelManager(void * stuff)=NULL);
     1746  /// Assignment operator. This copies the data
     1747  CoinPthreadStuff & operator=(const CoinPthreadStuff & rhs);
     1748  /// Destructor
     1749  ~CoinPthreadStuff (  );
     1750  /// set stop start
     1751  inline void setStopStart(int value)
     1752  { stopStart_=value;}
     1753#ifndef NUMBER_THREADS
     1754#define NUMBER_THREADS 8
    17121755#endif
     1756  // For waking up thread
     1757  inline pthread_mutex_t * mutexPointer(int which,int thread=0)
     1758  { return mutex_+which+3*thread;}
     1759  inline pthread_barrier_t * barrierPointer()
     1760  { return &barrier_;}
     1761  inline int whichLocked(int thread=0) const
     1762  { return locked_[thread];}
     1763  inline CoinThreadInfo * threadInfoPointer(int thread=0)
     1764  { return threadInfo_+thread;}
     1765  void startParallelTask(int type,int iThread,void * info=NULL);
     1766  int waitParallelTask(int type, int & iThread,bool allowIdle);
     1767  void waitAllTasks();
     1768  /// so thread can find out which one it is
     1769  int whichThread() const;
     1770  void sayIdle(int iThread);
     1771  //void startThreads(int numberThreads);
     1772  //void stopThreads();
     1773  // For waking up thread
     1774  pthread_mutex_t mutex_[3*(NUMBER_THREADS+1)];
     1775  pthread_barrier_t barrier_;
     1776  CoinThreadInfo threadInfo_[NUMBER_THREADS+1];
     1777  pthread_t abcThread_[NUMBER_THREADS+1];
     1778  int locked_[NUMBER_THREADS+1];
     1779  int stopStart_;
     1780  int numberThreads_;
     1781};
     1782void * clp_parallelManager(void * stuff);
     1783#endif
     1784#endif
  • trunk/Clp/src/ClpSimplexDual.cpp

    r2063 r2070  
    36313631               marker[1-iFlip][1] = CoinMin(marker[1-iFlip][1], numberPossiblySwapped);
    36323632
    3633                if (totalThru + thruThis >= fabs(dualOut_) ||
     3633               double check = fabs(totalThru+thruThis);
     3634               // add a bit
     3635               check += 1.0e-8+1.0e-10*check;
     3636               if (check >= fabs(dualOut_) ||
    36343637                         increaseInObjective + increaseInThis < 0.0) {
    36353638                    // We should be pivoting in this batch
     
    41664169     return status;
    41674170}
     4171//static int count_status=0;
     4172//static double obj_status=0.0;
     4173//static int check_status=123456789;//41754;
    41684174//static int count_alpha=0;
    41694175/* Checks if finished.  Updates status */
     
    43034309          //if ((count_alpha%5000)==0)
    43044310          //printf("count alpha %d\n",count_alpha);
     4311     }
     4312     if(progress_.infeasibility_[0]<1.0e-1 &&
     4313        primalTolerance_==1.0e-7&&progress_.iterationNumber_[0]>0&&
     4314        progress_.iterationNumber_[CLP_PROGRESS-1]-progress_.iterationNumber_[0]>25) {
     4315       // default - so user did not set
     4316       int iP;
     4317       double minAverage=COIN_DBL_MAX;
     4318       double maxAverage=0.0;
     4319       for (iP=0;iP<CLP_PROGRESS;iP++) {
     4320         int n=progress_.numberInfeasibilities_[iP];
     4321         if (!n) {
     4322           break;
     4323         } else {
     4324           double average=progress_.infeasibility_[iP];
     4325           if (average>0.1)
     4326             break;
     4327           average /= static_cast<double>(n);
     4328           minAverage=CoinMin(minAverage,average);
     4329           maxAverage=CoinMax(maxAverage,average);
     4330         }
     4331       }
     4332       if (iP==CLP_PROGRESS&&minAverage<1.0e-5&&maxAverage<1.0e-3) {
     4333         // change tolerance
     4334#if CBC_USEFUL_PRINTING>0
     4335         printf("CCchanging tolerance\n");
     4336#endif
     4337         primalTolerance_=1.0e-6;
     4338         dblParam_[ClpPrimalTolerance]=1.0e-6;
     4339         moreSpecialOptions_ |= 4194304;
     4340       }
    43054341     }
    43064342     // at this stage status is -3 or -4 if looks infeasible
     
    46914727     }
    46924728#if 0
     4729     count_status++;
     4730     if (!numberIterations_)
     4731       obj_status=-1.0e30;
     4732     if (objectiveValue()<obj_status-0.01) {
     4733       printf("Backward obj at %d from %g to %g\n",
     4734              count_status,obj_status,objectiveValue());
     4735     }
     4736     obj_status=objectiveValue();
     4737     if (count_status>=check_status-1) {
     4738       printf("Trouble ahead - count_status %d\n",count_status);
     4739     }
     4740#endif
     4741#if 0
    46934742     printf("IT %d %g %g(%d) %g(%d)\n",
    46944743            numberIterations_, objectiveValue(),
     
    47284777     /* If we are primal feasible and any dual infeasibilities are on
    47294778        free variables then it is better to go to primal */
    4730      if (!numberPrimalInfeasibilities_ && !numberDualInfeasibilitiesWithoutFree_ &&
    4731                numberDualInfeasibilities_)
     4779     if (!numberPrimalInfeasibilities_ && ((!numberDualInfeasibilitiesWithoutFree_ &&
     4780                                            numberDualInfeasibilities_)||
     4781                                           (moreSpecialOptions_&2097152)!=0))
    47324782          problemStatus_ = 10;
    47334783     // dual bound coming in
  • trunk/Clp/src/ClpSimplexOther.cpp

    r2028 r2070  
    12811281     return returnCode;
    12821282}
     1283/* Sets solution in dualized problem
     1284   non-zero return code indicates minor problems
     1285*/
     1286int
     1287ClpSimplexOther::setInDual(ClpSimplex * dualProblem)
     1288{
     1289  // Number of rows in dual problem was original number of columns
     1290  assert (numberColumns_ == dualProblem->numberRows());
     1291  // out If slack on d-row basic then column at bound otherwise column basic
     1292  // out If d-column basic then rhs tight
     1293  // if column at bound then slack on d-row basic
     1294  // if column basic then slack on d-row at bound
     1295  // if rhs non-basic then d-column basic
     1296  // if rhs basic then d-column ?
     1297  int numberBasic = 0;
     1298  int iRow, iColumn = 0;
     1299  //int numberExtraRows = dualProblem->numberColumns()-numberRows_;
     1300  //const double * objective = this->objective();
     1301  //double * dualDual = dualProblem->dualRowSolution();
     1302  //double * dualDj = dualProblem->dualColumnSolution();
     1303  double * dualSol = dualProblem->primalColumnSolution();
     1304  //double * dualActs = dualProblem->primalRowSolution();
     1305  const double * lower = dualProblem->columnLower();
     1306  const double * upper = dualProblem->columnUpper();
     1307  // position at bound information
     1308  int jColumn = numberRows_;
     1309  for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
     1310    Status status = getColumnStatus(iColumn);
     1311    Status statusD = dualProblem->getRowStatus(iColumn);
     1312    Status statusDJ = dualProblem->getColumnStatus(jColumn);
     1313    if (status==atLowerBound||
     1314        status==isFixed||
     1315        status==atUpperBound) {
     1316      dualProblem->setRowStatus(iColumn,basic);
     1317      numberBasic++;
     1318      if (columnUpper_[iColumn] < 1.0e20 &&
     1319          columnLower_[iColumn] > -1.0e20) {
     1320        bool mainLower =(fabs(columnLower_[iColumn]) < fabs(columnUpper_[iColumn]));
     1321        // fix this
     1322        if (mainLower) {
     1323          if (status==atUpperBound) {
     1324            dualProblem->setColumnStatus(jColumn,atUpperBound);
     1325          } else {
     1326            dualProblem->setColumnStatus(jColumn,atUpperBound);
     1327          }
     1328        } else {
     1329          if (status==atUpperBound) {
     1330            dualProblem->setColumnStatus(jColumn,atLowerBound);
     1331          } else {
     1332            dualProblem->setColumnStatus(jColumn,atLowerBound);
     1333          }
     1334        }
     1335        assert(statusDJ == dualProblem->getColumnStatus(jColumn));
     1336        jColumn++;
     1337      }
     1338    } else if (status==isFree) {
     1339      dualProblem->setRowStatus(iColumn,basic);
     1340      numberBasic++;
     1341    } else {
     1342      assert (status==basic);
     1343      //numberBasic++;
     1344    }
     1345    assert(statusD == dualProblem->getRowStatus(iColumn));
     1346  }
     1347  // now rows (no ranges at first)
     1348  for (iRow = 0; iRow < numberRows_; iRow++) {
     1349    Status status = getRowStatus(iRow);
     1350    Status statusD = dualProblem->getColumnStatus(iRow);
     1351    if (status == basic) {
     1352      // dual variable is at bound
     1353      if (!lower[iRow]) {
     1354        dualProblem->setColumnStatus(iRow,atLowerBound);
     1355      } else if (!upper[iRow]) {
     1356        dualProblem->setColumnStatus(iRow,atUpperBound);
     1357      } else {
     1358        dualProblem->setColumnStatus(iRow,isFree);
     1359        dualSol[iRow]=0.0;
     1360      }
     1361    } else {
     1362      // dual variable is basic
     1363      dualProblem->setColumnStatus(iRow,basic);
     1364      numberBasic++;
     1365    }
     1366    if (rowLower_[iRow] < -1.0e20 &&
     1367        rowUpper_[iRow] > 1.0e20) {
     1368      if (rowUpper_[iRow] != rowLower_[iRow]) {
     1369        printf("can't handle ranges yet\n");
     1370        abort();
     1371      }
     1372    }
     1373    assert(statusD == dualProblem->getColumnStatus(iRow));
     1374  }
     1375  if (numberBasic != numberColumns_) {
     1376    printf("Bad basis - ranges - coding needed ??\n");
     1377    abort();
     1378  }
     1379  return 0;
     1380}
    12831381/* Does very cursory presolve.
    12841382   rhs is numberRows, whichRows is 3*numberRows and whichColumns is 2*numberColumns
     
    15591657                         up[iRow] = CoinMax(upper - rowLower2[iRow], 0.0) + tolerance;
    15601658                    }
     1659                    // tighten row bounds
     1660                    if (lower>-1.0e10)
     1661                      rowLower2[iRow] = CoinMax(rowLower2[iRow],
     1662                                                lower - 1.0e-6*(1.0+fabs(lower)));
     1663                    if (upper<1.0e10)
     1664                      rowUpper2[iRow] = CoinMin(rowUpper2[iRow],
     1665                                                upper + 1.0e-6*(1.0+fabs(upper)));
    15611666               }
    15621667               if (!feasible) {
     
    28722977  ClpDataSave data = saveData();
    28732978  int numberTotal = numberRows_ + numberColumns_;
    2874   int ratio = (2*sizeof(int))/sizeof(double);
     2979  int ratio = static_cast<int>((2*sizeof(int))/sizeof(double));
    28752980  assert (ratio==1||ratio==2);
    28762981  // allow for unscaled - even if not needed
  • trunk/Clp/src/ClpSimplexOther.hpp

    r1884 r2070  
    226226         non-zero return code indicates minor problems
    227227     */
    228   int restoreFromDual(const ClpSimplex * dualProblem,
    229                       bool checkAccuracy=false);
     228     int restoreFromDual(const ClpSimplex * dualProblem,
     229                         bool checkAccuracy=false);
     230     /** Sets solution in dualized problem
     231         non-zero return code indicates minor problems
     232     */
     233     int setInDual(ClpSimplex * dualProblem);
    230234     /** Does very cursory presolve.
    231235         rhs is numberRows, whichRows is 3*numberRows and whichColumns is 2*numberColumns.
  • trunk/Clp/src/ClpSimplexPrimal.cpp

    r2030 r2070  
    908908                    problemStatus_ = -3;
    909909          }
     910          if(progress->realInfeasibility_[0]<1.0e-1 &&
     911             primalTolerance_==1.0e-7&&progress->iterationNumber_[0]>0&&
     912             progress->iterationNumber_[CLP_PROGRESS-1]-progress->iterationNumber_[0]>25) {
     913            // default - so user did not set
     914            int iP;
     915            double minAverage=COIN_DBL_MAX;
     916            double maxAverage=0.0;
     917            for (iP=0;iP<CLP_PROGRESS;iP++) {
     918              int n=progress->numberInfeasibilities_[iP];
     919              if (!n) {
     920                break;
     921              } else {
     922                double average=progress->realInfeasibility_[iP];
     923                if (average>0.1)
     924                  break;
     925                average /= static_cast<double>(n);
     926                minAverage=CoinMin(minAverage,average);
     927                maxAverage=CoinMax(maxAverage,average);
     928              }
     929            }
     930            if (iP==CLP_PROGRESS&&minAverage<1.0e-5&&maxAverage<1.0e-3) {
     931              // change tolerance
     932#if CBC_USEFUL_PRINTING>0
     933              printf("CCchanging tolerance\n");
     934#endif
     935              primalTolerance_=1.0e-6;
     936              dblParam_[ClpPrimalTolerance]=1.0e-6;
     937              moreSpecialOptions_ |= 4194304;
     938            }
     939          }
    910940          // at this stage status is -3 or -5 if looks unbounded
    911941          // get primal and dual solutions
  • trunk/Clp/src/ClpSolve.cpp

    r2044 r2070  
    777777      this->checkSolutionInternal();
    778778      if (sumDualInfeasibilities_>100.0*dualTolerance_) {
     779#if CBC_USEFUL_PRINTING>0
    779780        printf("internal check on duals failed %d %g\n",
    780781               numberDualInfeasibilities_,sumDualInfeasibilities_);
     782#endif
    781783      } else {
    782784        sumDualInfeasibilities_=0.0;
     
    784786      }
    785787      if (sumPrimalInfeasibilities_>100.0*primalTolerance_) {
     788#if CBC_USEFUL_PRINTING>0
    786789        printf("internal check on primals failed %d %g\n",
    787790               numberPrimalInfeasibilities_,sumPrimalInfeasibilities_);
     791#endif
    788792      } else {
    789793        sumPrimalInfeasibilities_=0.0;
     
    972976            //ClpModel::stopPermanentArrays();
    973977            //setSpecialOptions(specialOptions()&~65536);
     978            // try setting tolerances up
     979#define CLP_NEW_TOLERANCE 1.0e-7
     980            if (model2->primalTolerance()==1.0e-7&&model2->dualTolerance()==1.0e-7) {
     981              model2->setPrimalTolerance(CLP_NEW_TOLERANCE);
     982              model2->setDualTolerance(CLP_NEW_TOLERANCE);
     983            }
    974984#endif
    975985              model2->eventHandler()->setSimplex(model2);
     
    52255235  return model;
    52265236}
    5227 #ifdef ABC_INHERIT
    5228 // Use pthreads
    5229 #include <pthread.h>
    5230 class AbcPthreadStuff {
    5231 public:
    5232   /**@name Constructors and destructor and copy */
    5233   //@{
    5234   /** Main constructor
    5235   */
    5236   AbcPthreadStuff (int numberThreads=0);
    5237   /// Assignment operator. This copies the data
    5238   AbcPthreadStuff & operator=(const AbcPthreadStuff & rhs);
    5239   /// Destructor
    5240   ~AbcPthreadStuff (  );
    5241   /// set stop start
    5242   inline void setStopStart(int value)
    5243   { stopStart_=value;}
    5244 #ifndef NUMBER_THREADS
    5245 #define NUMBER_THREADS 8
    5246 #endif
    5247   // For waking up thread
    5248   inline pthread_mutex_t * mutexPointer(int which,int thread=0)
    5249   { return mutex_+which+3*thread;}
    5250   inline pthread_barrier_t * barrierPointer()
    5251   { return &barrier_;}
    5252   inline int whichLocked(int thread=0) const
    5253   { return locked_[thread];}
    5254   inline CoinAbcThreadInfo * threadInfoPointer(int thread=0)
    5255   { return threadInfo_+thread;}
    5256   void startParallelTask(int type,int iThread,void * info=NULL);
    5257   int waitParallelTask(int type, int & iThread,bool allowIdle);
    5258   void waitAllTasks();
    5259   /// so thread can find out which one it is
    5260   int whichThread() const;
    5261   //void startThreads(int numberThreads);
    5262   //void stopThreads();
    5263   // For waking up thread
    5264   pthread_mutex_t mutex_[3*(NUMBER_THREADS+1)];
    5265   pthread_barrier_t barrier_;
    5266   CoinAbcThreadInfo threadInfo_[NUMBER_THREADS+1];
    5267   pthread_t abcThread_[NUMBER_THREADS+1];
    5268   int locked_[NUMBER_THREADS+1];
    5269   int stopStart_;
    5270   int numberThreads_;
    5271 };
    5272 static void * abc_parallelManager(void * stuff)
     5237#if defined(ABC_INHERIT) || defined(CBC_THREAD)
     5238void * clp_parallelManager(void * stuff)
    52735239{
    5274   AbcPthreadStuff * driver = reinterpret_cast<AbcPthreadStuff *>(stuff);
     5240  CoinPthreadStuff * driver = reinterpret_cast<CoinPthreadStuff *>(stuff);
    52755241  int whichThread=driver->whichThread();
    5276   CoinAbcThreadInfo * threadInfo = driver->threadInfoPointer(whichThread);
     5242  CoinThreadInfo * threadInfo = driver->threadInfoPointer(whichThread);
    52775243  threadInfo->status=-1;
    52785244  int * which = threadInfo->stuff;
     
    53335299  }
    53345300}
    5335 AbcPthreadStuff::AbcPthreadStuff(int numberThreads)
     5301CoinPthreadStuff::CoinPthreadStuff(int numberThreads,
     5302                                   void * parallelManager(void * stuff))
    53365303{
    53375304  numberThreads_=numberThreads;
     
    53515318  pthread_barrier_init(&barrier_, /*&attr*/ NULL, numberThreads+1);
    53525319  for (int iThread=0;iThread<numberThreads;iThread++) {
    5353     pthread_create(&abcThread_[iThread], NULL, abc_parallelManager, reinterpret_cast<void *>(this));
     5320    pthread_create(&abcThread_[iThread], NULL, parallelManager, reinterpret_cast<void *>(this));
    53545321  }
    53555322  pthread_barrier_wait(&barrier_);
     
    53615328  }
    53625329}
    5363 AbcPthreadStuff::~AbcPthreadStuff()
     5330CoinPthreadStuff::~CoinPthreadStuff()
    53645331{
    53655332  for (int iThread=0;iThread<numberThreads_;iThread++) {
     
    53755342// so thread can find out which one it is
    53765343int
    5377 AbcPthreadStuff::whichThread() const
     5344CoinPthreadStuff::whichThread() const
    53785345{
    53795346  pthread_t thisThread=pthread_self();
     
    53875354}
    53885355void
    5389 AbcPthreadStuff::startParallelTask(int type, int iThread, void * info)
     5356CoinPthreadStuff::startParallelTask(int type, int iThread, void * info)
    53905357{
    53915358  /*
     
    54015368  pthread_mutex_unlock (&mutex_[locked_[iThread]+3*iThread]);
    54025369}
     5370void
     5371CoinPthreadStuff::sayIdle(int iThread)
     5372{
     5373  threadInfo_[iThread].status = -1;
     5374  threadInfo_[iThread].stuff[3] = -1;
     5375}
    54035376int
    5404 AbcPthreadStuff::waitParallelTask(int type,int & iThread, bool allowIdle)
     5377CoinPthreadStuff::waitParallelTask(int type,int & iThread, bool allowIdle)
    54055378{
    54065379  bool finished = false;
     
    54465419}
    54475420void
    5448 AbcPthreadStuff::waitAllTasks()
     5421CoinPthreadStuff::waitAllTasks()
    54495422{
    54505423  int nWait=0;
     
    58895862     //setAbcState(1);
    58905863     int numberCpu=CoinMin((this->abcState()&15),4);
    5891      AbcPthreadStuff threadInfo(numberCpu);
     5864     CoinPthreadStuff threadInfo(numberCpu,clp_parallelManager);
    58925865     masterModel.setAbcState(this->abcState());
    58935866     //AbcSimplex * tempMaster=masterModel.dealWithAbc(2,10,true);
  • trunk/Clp/src/ClpSolver.cpp

    r2030 r2070  
    844844#else
    845845                                   AbcSimplex * model2 = models + iModel;
     846#endif
    846847                                   if (type==CLP_PARAM_ACTION_EITHERSIMPLEX||
    847848                                       type==CBC_PARAM_ACTION_BAB)
    848849                                     solveOptions.setSpecialOption(3,0); // allow +-1
    849 #endif
    850850                                   if (dualize==4) {
    851851                                     solveOptions.setSpecialOption(4, 77);
  • trunk/Clp/src/Idiot.cpp

    r2030 r2070  
    615615     memset(lambda, 0, nrows * sizeof(double));
    616616     slackStart = countCostedSlacks(model_);
     617     // redo in case odd matrix
     618     row = matrix->getIndices();
     619     columnStart = matrix->getVectorStarts();
     620     columnLength = matrix->getVectorLengths();
     621     element = matrix->getElements();
    617622     if (slackStart >= 0) {
    618623       COIN_DETAIL_PRINT(printf("This model has costed slacks\n"));
  • trunk/Clp/src/OsiClp/OsiClpSolverInterface.cpp

    r2063 r2070  
    3636#include "ClpPresolve.hpp"
    3737#include "CoinLpIO.hpp"
     38//#define PRINT_TIME
     39#ifdef PRINT_TIME
    3840static double totalTime=0.0;
     41#endif
    3942//#define SAVE_MODEL 1
    4043#ifdef SAVE_MODEL
     
    7275  bool deleteSolver;
    7376  ClpSimplex * solver;
     77#ifdef PRINT_TIME
    7478  double time1 = CoinCpuTime();
     79#endif
    7580  int userFactorizationFrequency = modelPtr_->factorization()->maximumPivots();
    7681  int totalIterations=0;
     
    209214    // sort out hints;
    210215    // algorithm 0 whatever, -1 force dual, +1 force primal
    211     int algorithm = 0;
     216    int algorithm = 0; 
    212217    gotHint = (getHintParam(OsiDoDualInInitial,takeHint,strength));
    213218    assert (gotHint);
     
    769774  // mark so we can pick up objective value quickly
    770775  modelPtr_->upperIn_=0.0;
     776#ifdef PRINT_TIME
    771777  time1 = CoinCpuTime()-time1;
    772778  totalTime += time1;
     779#endif
    773780  assert (!modelPtr_->disasterHandler());
    774781  if (lastAlgorithm_<1||lastAlgorithm_>2)
     
    786793  modelPtr_->rowCopy_=NULL;
    787794#endif
    788   //std::cout<<time1<<" seconds - total "<<totalTime<<std::endl;
     795#ifdef PRINT_TIME
     796  std::cout<<time1<<" seconds - total "<<totalTime<<std::endl;
     797#endif
    789798  delete pinfo;
    790799}
     
    10911100      }
    10921101#endif
    1093       if((specialOptions_&1)==0||(specialOptions_&2048)!=0/*||skipCrunch*/) {
     1102      if((specialOptions_&1)==0||
     1103         (specialOptions_&2048)!=0||
     1104         (modelPtr_->specialOptions_&2097152)!=0/*||skipCrunch*/) {
    10941105        disasterHandler_->setWhereFrom(0); // dual
    10951106        if (inCbcOrOther)
     
    14181429    break;
    14191430  }
     1431  //#define NO_CRUNCH2
     1432#ifdef NO_CRUNCH2
     1433  specialOptions_ &= ~1;
     1434#endif
    14201435  bool stopPrinting=false;
    14211436  if (printOut<0) {
     
    22392254      }
    22402255    }
    2241     //#define REDO_STEEPEST_EDGE
    2242 #ifdef REDO_STEEPEST_EDGE
    2243     if (dynamic_cast<ClpDualRowSteepest *>(modelPtr_->dualRowPivot())) {
    2244       // new copy (should think about keeping around)
    2245       ClpDualRowSteepest steepest;
    2246       modelPtr_->setDualRowPivotAlgorithm(steepest);
    2247     }     
    2248 #endif
    22492256    // Start of fast iterations
    22502257    bool alwaysFinish= ((specialOptions_&32)==0) ? true : false;
     
    22622269    if (problemStatus==10||problemStatus<0) {
    22632270      // was trying to clean up or something odd
    2264       if (problemStatus==10)
     2271      if (problemStatus==10) {
     2272        obj = saveObjectiveValue;
    22652273        lastAlgorithm_=1; // so won't fail on cutoff (in CbcNode)
     2274      }
    22662275      status=1;
    22672276    }
     
    22872296        } else if (problemStatus==10) {
    22882297          problemStatus=3;
     2298          obj = saveObjectiveValue;
    22892299        } else if (!modelPtr_->numberPrimalInfeasibilities()) {
    22902300          problemStatus=1; // infeasible
     
    22952305        //modelPtr_->computeObjectiveValue();
    22962306        lastAlgorithm_=1; // so won't fail on cutoff (in CbcNode)
     2307        obj = saveObjectiveValue;
    22972308        problemStatus=3;
    22982309      }
     
    24102421    double * lowerSmall = smallModel_->lowerRegion();
    24112422    double * upperSmall = smallModel_->upperRegion();
     2423    double * solutionSmall = smallModel_->solutionRegion();
    24122424    double * lowerSmallReal = smallModel_->columnLower();
    24132425    double * upperSmallReal = smallModel_->columnUpper();
     
    24292441          value /= columnScale[i];
    24302442        lowerSmall[i]=value;
     2443        if (smallModel_->getColumnStatus(i)==ClpSimplex::atLowerBound)
     2444          solutionSmall[i]=value;
    24312445      }
    24322446      if (upperBig[iColumn]<saveUpperOriginal[iColumn]) {
     
    24372451          value /= columnScale[i];
    24382452        upperSmall[i]=value;
     2453        if (smallModel_->getColumnStatus(i)==ClpSimplex::atUpperBound)
     2454          solutionSmall[i]=value;
    24392455      }
    24402456      if (upperSmall[i]<lowerSmall[i]-1.0e-8)
     
    24522468    int status;
    24532469    if (!infeasible) {
    2454 #ifdef REDO_STEEPEST_EDGE
    2455       if (dynamic_cast<ClpDualRowSteepest *>(smallModel_->dualRowPivot())) {
    2456         // new copy
    2457         ClpDualRowSteepest steepest;
    2458         smallModel_->setDualRowPivotAlgorithm(steepest);
    2459       }     
    2460 #endif
     2470      if((modelPtr_->specialOptions()&0x011200000)==0x11200000) {
     2471        smallModel_->specialOptions_|=2097152;
     2472        //smallModel_->specialOptions_&=~2097152;
     2473        delete [] modelPtr_->ray_;
     2474        delete [] smallModel_->ray_;
     2475        modelPtr_->ray_=NULL;
     2476        smallModel_->ray_=NULL;
     2477      }
    24612478      status = static_cast<ClpSimplexDual *>(smallModel_)->fastDual(alwaysFinish);
     2479      if((modelPtr_->specialOptions()&0x011200000)==0x11200000&&smallModel_->ray_) {
     2480        if (smallModel_->sequenceOut_<smallModel_->numberColumns_)
     2481            modelPtr_->sequenceOut_ = whichColumn[smallModel_->sequenceOut_];
     2482        else
     2483            modelPtr_->sequenceOut_ = whichRow[smallModel_->sequenceOut_-smallModel_->numberColumns_]+modelPtr_->numberColumns_;
     2484        if (smallModel_->rowScale_) {
     2485          // scale ray
     2486          double scaleFactor=1.0;
     2487          if (smallModel_->sequenceOut_<smallModel_->numberColumns_)
     2488            scaleFactor = smallModel_->columnScale_[smallModel_->sequenceOut_];
     2489          for (int i = 0; i < smallModel_->numberRows_; i++) {
     2490            smallModel_->ray_[i] *= smallModel_->rowScale_[i]*scaleFactor;
     2491          }
     2492        }
     2493      }
    24622494    } else {
    24632495      status=0;
     
    25032535    if (problemStatus==10||problemStatus<0) {
    25042536      // was trying to clean up or something odd
    2505       if (problemStatus==10)
     2537      if (problemStatus==10) 
    25062538        lastAlgorithm_=1; // so won't fail on cutoff (in CbcNode)
    25072539      status=1;
     
    25962628    if (modelPtr_==smallModel_)
    25972629      modelPtr_->computeObjectiveValue();
     2630    if (problemStatus==1&&smallModel_->ray_) {
     2631      delete [] modelPtr_->ray_;
     2632      // get ray to full problem
     2633      int numberRows = modelPtr_->numberRows();
     2634      int numberRows2 = smallModel_->numberRows();
     2635      int nCopy = 3*numberRows+2*numberColumns;
     2636      int nBound = whichRow[nCopy];
     2637      double * ray = new double [numberRows];
     2638      memset(ray,0,numberRows*sizeof(double));
     2639      for (int i = 0; i < numberRows2; i++) {
     2640        int iRow = whichRow[i];
     2641        ray[iRow] = smallModel_->ray_[i];
     2642      }
     2643      // Column copy of matrix
     2644      const double * element = getMatrixByCol()->getElements();
     2645      const int * row = getMatrixByCol()->getIndices();
     2646      const CoinBigIndex * columnStart = getMatrixByCol()->getVectorStarts();
     2647      const int * columnLength = getMatrixByCol()->getVectorLengths();
     2648      // translate
     2649      for (int jRow = nBound; jRow < 2 * numberRows; jRow++) {
     2650        int iRow = whichRow[jRow];
     2651        int iColumn = whichRow[jRow+numberRows];
     2652        if (modelPtr_->getColumnStatus(iColumn) == ClpSimplex::basic) {
     2653          double value = 0.0;
     2654          double sum = 0.0;
     2655          for (CoinBigIndex j = columnStart[iColumn];
     2656               j < columnStart[iColumn] + columnLength[iColumn]; j++) {
     2657            if (iRow == row[j]) {
     2658              value = element[j];
     2659            } else {
     2660              sum += ray[row[j]]*element[j];
     2661            }
     2662          }
     2663          ray[iRow] = -sum / value;
     2664        }
     2665      }
     2666      for (int i=0;i<modelPtr_->numberColumns_;i++) {
     2667        if (modelPtr_->getStatus(i)!=ClpSimplex::basic&&
     2668            modelPtr_->columnLower_[i]==modelPtr_->columnUpper_[i])
     2669          modelPtr_->setStatus(i,ClpSimplex::isFixed);
     2670      }
     2671      modelPtr_->ray_=ray;
     2672      //delete [] smallModel_->ray_;
     2673      //smallModel_->ray_=NULL;
     2674      modelPtr_->directionOut_=smallModel_->directionOut_;
     2675      lastAlgorithm_=2; // dual
     2676    }
    25982677#if 1
    25992678    if (status&&!problemStatus) {
     
    27052784}
    27062785
     2786#ifdef CONFLICT_CUTS
     2787// Return a conflict analysis cut from small model
     2788OsiRowCut *
     2789OsiClpSolverInterface::smallModelCut(const double * originalLower, const double * originalUpper,
     2790                                     int numberRowsAtContinuous,const int * whichGenerator,
     2791                                     int typeCut)
     2792{
     2793  if(smallModel_&&smallModel_->ray_) {
     2794    //printf("Could do small cut\n");
     2795    int numberRows = modelPtr_->numberRows();
     2796    int numberRows2 = smallModel_->numberRows();
     2797    int numberColumns = modelPtr_->numberColumns();
     2798    int numberColumns2 = smallModel_->numberColumns();
     2799    double * arrayD = reinterpret_cast<double *> (spareArrays_);
     2800    double * saveSolution = arrayD+1;
     2801    double * saveLower = saveSolution + (numberRows+numberColumns);
     2802    double * saveUpper = saveLower + (numberRows+numberColumns);
     2803    double * saveObjective = saveUpper + (numberRows+numberColumns);
     2804    double * saveLowerOriginal = saveObjective + (numberRows+numberColumns);
     2805    double * saveUpperOriginal = saveLowerOriginal + numberColumns;
     2806    //double * lowerOriginal = modelPtr_->columnLower();
     2807    //double * upperOriginal = modelPtr_->columnUpper();
     2808    int * savePivot = reinterpret_cast<int *> (saveUpperOriginal + numberColumns);
     2809    int * whichRow = savePivot+numberRows;
     2810    int * whichColumn = whichRow + 3*numberRows;
     2811    int nCopy = 3*numberRows+2*numberColumns;
     2812    int nBound = whichRow[nCopy];
     2813    if (smallModel_->sequenceOut_>=0&&smallModel_->sequenceOut_<numberColumns2)
     2814      modelPtr_->sequenceOut_ = whichColumn[smallModel_->sequenceOut_];
     2815    else
     2816      modelPtr_->sequenceOut_ = modelPtr_->numberColumns_+whichRow[smallModel_->sequenceOut_];
     2817    unsigned char * saveStatus = CoinCopyOfArray(modelPtr_->status_,
     2818                                                 numberRows+numberColumns);
     2819    // get ray to full problem
     2820    for (int i = 0; i < numberColumns2; i++) {
     2821      int iColumn = whichColumn[i];
     2822      modelPtr_->setStatus(iColumn, smallModel_->getStatus(i));
     2823    }
     2824    double * ray = new double [numberRows+numberColumns2+numberColumns];
     2825    char * mark = new char [numberRows];
     2826    memset(ray,0,(numberRows+numberColumns2+numberColumns)*sizeof(double));
     2827    double * smallFarkas = ray+numberRows;
     2828    double * farkas = smallFarkas+numberColumns2;
     2829    double * saveScale = smallModel_->rowScale_;
     2830    smallModel_->rowScale_=NULL;
     2831    smallModel_->transposeTimes(1.0,smallModel_->ray_,smallFarkas);
     2832    smallModel_->rowScale_=saveScale;
     2833    for (int i=0;i<numberColumns2;i++)
     2834      farkas[whichColumn[i]]=smallFarkas[i];
     2835    memset(mark,0,numberRows);
     2836    for (int i = 0; i < numberRows2; i++) {
     2837      int iRow = whichRow[i];
     2838      modelPtr_->setRowStatus(iRow, smallModel_->getRowStatus(i));
     2839      ray[iRow] = smallModel_->ray_[i];
     2840      mark[iRow]=1;
     2841    }
     2842    // Column copy of matrix
     2843    const double * element = getMatrixByCol()->getElements();
     2844    const int * row = getMatrixByCol()->getIndices();
     2845    const CoinBigIndex * columnStart = getMatrixByCol()->getVectorStarts();
     2846    const int * columnLength = getMatrixByCol()->getVectorLengths();
     2847    // pick up small pivot row
     2848    int pivotRow = smallModel_->spareIntArray_[3];
     2849    // translate
     2850    if (pivotRow>=0)
     2851      pivotRow=whichRow[pivotRow];
     2852    modelPtr_->spareIntArray_[3]=pivotRow;
     2853    for (int jRow = nBound; jRow < 2 * numberRows; jRow++) {
     2854      int iRow = whichRow[jRow];
     2855      int iColumn = whichRow[jRow+numberRows];
     2856      if (modelPtr_->getColumnStatus(iColumn) == ClpSimplex::basic) {
     2857        double value = 0.0;
     2858        double sum = 0.0;
     2859        for (CoinBigIndex j = columnStart[iColumn];
     2860             j < columnStart[iColumn] + columnLength[iColumn]; j++) {
     2861          if (iRow == row[j]) {
     2862            value = element[j];
     2863          } else if (mark[row[j]]) {
     2864            sum += ray[row[j]]*element[j];
     2865          }
     2866        }
     2867        double target=farkas[iColumn];
     2868        if (iRow!=pivotRow) {
     2869          ray[iRow] = (target-sum) / value;
     2870        } else {
     2871          printf("what now - direction %d wanted %g sum %g value %g\n",
     2872                 smallModel_->directionOut_,ray[iRow],
     2873                 sum,value);
     2874        }
     2875        mark[iRow]=1;
     2876      }
     2877    }
     2878    delete [] mark;
     2879    for (int i=0;i<modelPtr_->numberColumns_;i++) {
     2880      if (modelPtr_->getStatus(i)!=ClpSimplex::basic&&
     2881          modelPtr_->columnLower_[i]==modelPtr_->columnUpper_[i])
     2882        modelPtr_->setStatus(i,ClpSimplex::isFixed);
     2883    }
     2884#if 0
     2885    {
     2886      int nBad=0;
     2887      double * fBig = new double [2*numberColumns+2*numberRows];
     2888      double * fSmall = fBig+numberColumns;
     2889      double * effectiveRhs = fSmall+numberColumns;
     2890      double * effectiveRhs2 = effectiveRhs+numberRows;
     2891      memset(fBig,0,2*numberColumns*sizeof(double));
     2892      {
     2893        memcpy(fBig,modelPtr_->columnLower_,numberColumns*sizeof(double));
     2894        const double * columnLower = modelPtr_->columnLower_;
     2895        const double * columnUpper = modelPtr_->columnUpper_;
     2896        for (int i=0;i<numberColumns2;i++) {
     2897          int iBig=whichColumn[i];
     2898          if (smallModel_->getStatus(i)==ClpSimplex::atLowerBound||
     2899              smallModel_->getStatus(i)==ClpSimplex::isFixed||
     2900              columnLower[iBig]==columnUpper[iBig]) {
     2901            fSmall[i]=columnLower[iBig];
     2902          } else if (smallModel_->getStatus(i)==ClpSimplex::atUpperBound) {
     2903            fSmall[i]=columnUpper[iBig];
     2904          } else if (i==smallModel_->sequenceOut_) {
     2905            if (smallModel_->directionOut_<0) {
     2906              // above upper bound
     2907              fSmall[i]=columnUpper[iBig];
     2908            } else {
     2909              // below lower bound
     2910              fSmall[i]=columnLower[iBig];
     2911            }
     2912          } else if (smallModel_->columnLower_[i]==smallModel_->columnUpper_[i]) {
     2913            fSmall[i]=smallModel_->columnLower_[i];
     2914          }
     2915          fBig[whichColumn[i]]=fSmall[i];
     2916        }
     2917        {
     2918          // only works for one test case
     2919          memcpy(effectiveRhs2,modelPtr_->rowUpper_,numberRows*sizeof(double));
     2920          double * saveScale = modelPtr_->rowScale_;
     2921          ClpPackedMatrix * saveMatrix = modelPtr_->scaledMatrix_;
     2922          modelPtr_->rowScale_=NULL;
     2923          modelPtr_->scaledMatrix_=NULL;
     2924          modelPtr_->times(-1.0,fBig,effectiveRhs2);
     2925          modelPtr_->scaledMatrix_=saveMatrix;
     2926          modelPtr_->rowScale_=saveScale;
     2927          memset(fBig,0,numberColumns*sizeof(double));
     2928        }
     2929        const double * rowLower = smallModel_->rowLower_;
     2930        const double * rowUpper = smallModel_->rowUpper_;
     2931        int pivotRow = smallModel_->spareIntArray_[3];
     2932        for (int i=0;i<numberRows2;i++) {
     2933          if (smallModel_->getRowStatus(i)==ClpSimplex::atLowerBound||
     2934              rowUpper[i]==rowLower[i]||
     2935              smallModel_->getRowStatus(i)==ClpSimplex::isFixed) {
     2936            effectiveRhs[i]=rowLower[i];
     2937            //if (smallModel_->getRowStatus(i)==ClpSimplex::atLowerBound)
     2938            //assert (ray[i]>-1.0e-3);
     2939          } else if (smallModel_->getRowStatus(i)==ClpSimplex::atUpperBound) {
     2940            effectiveRhs[i]=rowUpper[i];
     2941            //assert (ray[i]<1.0e-3);
     2942          } else {
     2943            if (smallModel_->getRowStatus(i)!=ClpSimplex::basic) {
     2944              assert (rowUpper[i]>1.0e30||rowLower[i]<-1.0e30); // eventually skip
     2945              if (rowUpper[i]<1.0e30)
     2946                effectiveRhs[i]=rowUpper[i];
     2947              else
     2948                effectiveRhs[i]=rowLower[i];
     2949            }
     2950          }
     2951          if (smallModel_->getRowStatus(i)==ClpSimplex::basic) {
     2952            effectiveRhs[i]=0.0;
     2953            // we should leave pivot row in and use direction for bound
     2954            if (fabs(smallModel_->ray_[i])>1.0e-8) {
     2955              printf("Basic small slack value %g on %d - pivotRow %d\n",smallModel_->ray_[i],i,pivotRow);
     2956              if (i==pivotRow) {
     2957                if (smallModel_->directionOut_<0)
     2958                  effectiveRhs[i]=rowUpper[i];
     2959                else
     2960                  effectiveRhs[i]=rowLower[i];
     2961              } else {
     2962                //smallModel_->ray_[i]=0.0;
     2963              }
     2964            }
     2965          }
     2966        }
     2967        //ClpPackedMatrix * saveMatrix = smallModel_->scaledMatrix_;
     2968        double * saveScale = smallModel_->rowScale_;
     2969        smallModel_->rowScale_=NULL;
     2970        smallModel_->scaledMatrix_=NULL;
     2971        smallModel_->times(-1.0,fSmall,effectiveRhs);
     2972        double bSum=0.0;
     2973        for (int i=0;i<numberRows2;i++) {
     2974          bSum += effectiveRhs[i]*smallModel_->ray_[i];
     2975        }
     2976        printf("Small bsum %g\n",bSum);
     2977        smallModel_->rowScale_=saveScale;
     2978      }
     2979      double * saveScale = smallModel_->rowScale_;
     2980      smallModel_->rowScale_=NULL;
     2981      memset(fSmall,0,numberColumns*sizeof(double));
     2982      smallModel_->transposeTimes(1.0,smallModel_->ray_,fSmall);
     2983      smallModel_->rowScale_=saveScale;
     2984      ClpPackedMatrix * saveMatrix = modelPtr_->scaledMatrix_;
     2985      saveScale = modelPtr_->rowScale_;
     2986      modelPtr_->rowScale_=NULL;
     2987      modelPtr_->scaledMatrix_=NULL;
     2988      modelPtr_->transposeTimes(1.0,ray,fBig);
     2989      modelPtr_->scaledMatrix_=saveMatrix;
     2990      modelPtr_->rowScale_=saveScale;
     2991      for (int i=0;i<numberColumns2;i++) {
     2992        int iBig=whichColumn[i];
     2993        if (fabs(fBig[iBig]-fSmall[i])>1.0e-4) {
     2994          printf("small %d %g big %d %g\n",
     2995                 i,fSmall[i],iBig,fBig[iBig]);
     2996          nBad++;
     2997        }
     2998      }
     2999      delete [] fBig;
     3000      if (nBad)
     3001        printf("Bad %d odd %d\n",nBad,2*numberRows-nBound);
     3002    }
     3003#endif
     3004    modelPtr_->ray_=ray;
     3005    lastAlgorithm_=2;
     3006    modelPtr_->directionOut_=smallModel_->directionOut_;
     3007    OsiRowCut * cut = modelCut(originalLower,originalUpper,
     3008                               numberRowsAtContinuous,whichGenerator,typeCut);
     3009    smallModel_->deleteRay();
     3010    memcpy(modelPtr_->status_,saveStatus,numberRows+numberColumns);
     3011    delete [] saveStatus;
     3012    if (cut) {
     3013      //printf("got small cut\n");
     3014      //delete cut;
     3015      //cut=NULL;
     3016    }
     3017    return cut;
     3018  } else {
     3019    return NULL;
     3020  }
     3021}
     3022static int debugMode=0;
     3023// Return a conflict analysis cut from model
     3024//      If type is 0 then genuine cut, if 1 then only partially processed
     3025OsiRowCut *
     3026OsiClpSolverInterface::modelCut(const double * originalLower, const double * originalUpper,
     3027                                int numberRowsAtContinuous,const int * whichGenerator,
     3028                                int typeCut)
     3029{
     3030  OsiRowCut * cut=NULL;
     3031  //return NULL;
     3032  if(modelPtr_->ray_) {
     3033    if (lastAlgorithm_==2) {
     3034      //printf("Could do normal cut\n");
     3035      // could use existing arrays
     3036      int numberRows=modelPtr_->numberRows_;
     3037      int numberColumns=modelPtr_->numberColumns_;
     3038      double * farkas = new double [2*numberColumns+numberRows];
     3039      double * bound = farkas + numberColumns;
     3040      double * effectiveRhs = bound + numberColumns;
     3041      /*const*/ double * ray = modelPtr_->ray_;
     3042      // have to get rid of local cut rows
     3043      whichGenerator -= numberRowsAtContinuous;
     3044      int badRows=0;
     3045      for (int iRow=numberRowsAtContinuous;iRow<numberRows;iRow++) {
     3046        int iType=whichGenerator[iRow];
     3047        if ((iType>=0&&iType<20000)) {
     3048          if (fabs(ray[iRow])>1.0e-10) {
     3049            badRows++;
     3050          }
     3051          ray[iRow]=0.0;
     3052        }
     3053      }
     3054      ClpSimplex debugModel;
     3055      if ((debugMode&4)!=0)
     3056        debugModel = *modelPtr_;
     3057      if (badRows&&(debugMode&1)!=0)
     3058        printf("%d rows from local cuts\n",badRows);
     3059      // get farkas row
     3060      ClpPackedMatrix * saveMatrix = modelPtr_->scaledMatrix_;
     3061      double * saveScale = modelPtr_->rowScale_;
     3062      modelPtr_->rowScale_=NULL;
     3063      modelPtr_->scaledMatrix_=NULL;
     3064      memset(farkas,0,(2*numberColumns+numberRows)*sizeof(double));
     3065      modelPtr_->transposeTimes(-1.0,ray,farkas);
     3066      //const char * integerInformation = modelPtr_->integerType_;
     3067      //assert (integerInformation);
     3068
     3069      int sequenceOut = modelPtr_->sequenceOut_;
     3070      // Put nonzero bounds in bound
     3071      const double * columnLower = modelPtr_->columnLower_;
     3072      const double * columnUpper = modelPtr_->columnUpper_;
     3073      const double * columnValue = modelPtr_->columnActivity_;
     3074      int numberBad=0;
     3075      int nNonzeroBasic=0;
     3076      for (int i=0;i<numberColumns;i++) {
     3077        double value = farkas[i];
     3078        double boundValue=0.0;
     3079        if (modelPtr_->getStatus(i)==ClpSimplex::basic) {
     3080          // treat as zero if small
     3081          if (fabs(value)<1.0e-8) {
     3082            value=0.0;
     3083            farkas[i]=0.0;
     3084          }
     3085          if (value) {
     3086            //printf("basic %d direction %d farkas %g\n",
     3087            //i,modelPtr_->directionOut_,value);
     3088            nNonzeroBasic++;
     3089            if (value<0.0)
     3090              boundValue=columnLower[i];
     3091            else
     3092              boundValue=columnUpper[i];
     3093          }
     3094        } else if (fabs(value)>1.0e-10) {
     3095          if (value<0.0) {
     3096            if (columnValue[i]>columnLower[i]+1.0e-5&&value<-1.0e-7) {
     3097              //printf("bad %d alpha %g not at lower\n",i,value);
     3098              numberBad++;
     3099            }
     3100            boundValue=columnLower[i];
     3101          } else {
     3102            if (columnValue[i]<columnUpper[i]-1.0e-5&&value>1.0e-7) {
     3103              //printf("bad %d alpha %g not at upper\n",i,value);
     3104              numberBad++;
     3105            }
     3106            boundValue=columnUpper[i];
     3107          }
     3108        }
     3109        bound[i]=boundValue;
     3110        if (fabs(boundValue)>1.0e10)
     3111          numberBad++;
     3112      }
     3113      const double * rowLower = modelPtr_->rowLower_;
     3114      const double * rowUpper = modelPtr_->rowUpper_;
     3115      //int pivotRow = modelPtr_->spareIntArray_[3];
     3116      //bool badPivot=pivotRow<0;
     3117      for (int i=0;i<numberRows;i++) {
     3118        double value = ray[i];
     3119        double rhsValue=0.0;
     3120        if (modelPtr_->getRowStatus(i)==ClpSimplex::basic) {
     3121          // treat as zero if small
     3122          if (fabs(value)<1.0e-8) {
     3123            value=0.0;
     3124            ray[i]=0.0;
     3125          }
     3126          if (value) {
     3127            //printf("row basic %d direction %d ray %g\n",
     3128            //     i,modelPtr_->directionOut_,value);
     3129            nNonzeroBasic++;
     3130            if (value<0.0)
     3131              rhsValue=rowLower[i];
     3132            else
     3133              rhsValue=rowUpper[i];
     3134          }
     3135        } else if (fabs(value)>1.0e-10) {
     3136          if (value<0.0)
     3137            rhsValue=rowLower[i];
     3138          else
     3139            rhsValue=rowUpper[i];
     3140        }
     3141        effectiveRhs[i]=rhsValue;
     3142      }
     3143      {
     3144        double bSum=0.0;
     3145        for (int i=0;i<numberRows;i++) {
     3146          bSum += effectiveRhs[i]*ray[i];
     3147        }
     3148        //printf("before bounds - bSum %g\n",bSum);
     3149      }
     3150      modelPtr_->times(-1.0,bound,effectiveRhs);
     3151      double bSum=0.0;
     3152      for (int i=0;i<numberRows;i++) {
     3153        bSum += effectiveRhs[i]*ray[i];
     3154      }
     3155#if 0
     3156      double bSum2=0.0;
     3157#if 1
     3158      {
     3159        int iSequence = modelPtr_->sequenceOut_;
     3160        double lower,value,upper;
     3161        if (iSequence<numberColumns) {
     3162          lower=columnLower[iSequence];
     3163          value=columnValue[iSequence];
     3164          upper=columnUpper[iSequence];
     3165        } else {
     3166          iSequence -= numberColumns;
     3167          lower=rowLower[iSequence];
     3168          value=modelPtr_->rowActivity_[iSequence];
     3169          upper=rowUpper[iSequence];
     3170        }
     3171        if (value<lower)
     3172          bSum2=value-lower;
     3173        else if (value>upper)
     3174          bSum2=-(value-upper);
     3175        else
     3176          printf("OUCHXX %g %g %g\n",
     3177                 lower,value,upper);
     3178      }
     3179#else
     3180      if (modelPtr_->valueOut_<modelPtr_->lowerOut_)
     3181        bSum2=modelPtr_->valueOut_-modelPtr_->lowerOut_;
     3182      else if (modelPtr_->valueOut_>modelPtr_->upperOut_)
     3183        bSum2=-(modelPtr_->valueOut_-modelPtr_->upperOut_);
     3184      else
     3185        printf("OUCHXX %g %g %g\n",
     3186               modelPtr_->lowerOut_,modelPtr_->valueOut_,modelPtr_->upperOut_);
     3187#endif
     3188#if 0
     3189      if (fabs(bSum-bSum2)>1.0e-5)
     3190        printf("XX ");
     3191      printf("after bounds - bSum %g - bSum2 %g\n",bSum,bSum2);
     3192#endif
     3193#endif
     3194      modelPtr_->scaledMatrix_=saveMatrix;
     3195      modelPtr_->rowScale_=saveScale;
     3196      if (numberBad||bSum>-1.0e-4/*||nNonzeroBasic>1*/ /*||bSum2>-1.0e-4*/) {
     3197#if PRINT_CONFLICT>1 //ndef NDEBUG
     3198        printf("bad BOUND bSum %g (bSum2 %g) - %d bad - %d basic\n",
     3199               bSum,bSum2,numberBad,nNonzeroBasic);
     3200#endif
     3201      } else {
     3202        if (numberColumns<0)
     3203          debugMode=-numberColumns;
     3204        if ((debugMode&4)!=0) {
     3205          int * tempC = new int[numberColumns];
     3206          double * temp = new double[numberColumns];
     3207          int n=0;
     3208          for (int i=0;i<numberColumns;i++) {
     3209            if (fabs(farkas[i])>1.0e-12) {
     3210              temp[n]=farkas[i];
     3211              tempC[n++]=i;
     3212            }
     3213          }
     3214          debugModel.addRow(n,tempC,temp,bSum,bSum);
     3215          delete [] tempC;
     3216          delete [] temp;
     3217        }
     3218        if ((debugMode&2)!=0) {
     3219          ClpSimplex dual=*modelPtr_;
     3220          dual.setLogLevel(63);
     3221          dual.scaling(0);
     3222          dual.dual();
     3223          assert (dual.problemStatus_==1);
     3224          if (dual.ray_) {
     3225            double * farkas2 = dual.reducedCost_;
     3226            // Put nonzero bounds in farkas
     3227            const double * columnLower = dual.columnLower_;
     3228            const double * columnUpper = dual.columnUpper_;
     3229            for (int i=0;i<numberColumns;i++) {
     3230              farkas2[i]=0.0;
     3231              if (dual.getStatus(i)==ClpSimplex::atLowerBound||
     3232                  columnLower[i]==columnUpper[i]) {
     3233                farkas2[i]=columnLower[i];
     3234              } else if (dual.getStatus(i)==ClpSimplex::atUpperBound) {
     3235                farkas2[i]=columnUpper[i];
     3236              } else if (i==dual.sequenceOut_) {
     3237                if (dual.directionOut_<0) {
     3238                  // above upper bound
     3239                  farkas2[i]=columnUpper[i];
     3240                } else {
     3241                  // below lower bound
     3242                  farkas2[i]=columnLower[i];
     3243                }
     3244              } else if (columnLower[i]==columnUpper[i]) {
     3245                farkas2[i]=columnLower[i];
     3246              }
     3247            }
     3248            double * effectiveRhs2 = dual.rowActivity_;
     3249            const double * rowLower = dual.rowLower_;
     3250            const double * rowUpper = dual.rowUpper_;
     3251            int pivotRow = dual.spareIntArray_[3];
     3252            for (int i=0;i<numberRows;i++) {
     3253              if (dual.getRowStatus(i)==ClpSimplex::atLowerBound||
     3254                  rowUpper[i]==rowLower[i]||
     3255                  dual.getRowStatus(i)==ClpSimplex::isFixed) {
     3256                effectiveRhs2[i]=rowLower[i];
     3257              } else if (dual.getRowStatus(i)==ClpSimplex::atUpperBound) {
     3258                effectiveRhs2[i]=rowUpper[i];
     3259              } else {
     3260                if (dual.getRowStatus(i)!=ClpSimplex::basic) {
     3261                  assert (rowUpper[i]>1.0e30||rowLower[i]<-1.0e30); // eventually skip
     3262                  if (rowUpper[i]<1.0e30)
     3263                    effectiveRhs2[i]=rowUpper[i];
     3264                  else
     3265                    effectiveRhs2[i]=rowLower[i];
     3266                }
     3267              }
     3268              if (dual.getRowStatus(i)==ClpSimplex::basic) {
     3269                effectiveRhs2[i]=0.0;
     3270                // we should leave pivot row in and use direction for bound
     3271                if (fabs(dual.ray_[i])>1.0e-8) {
     3272                  printf("Basic slack value %g on %d - pivotRow %d\n",ray[i],i,pivotRow);
     3273                  if (i==pivotRow) {
     3274                    if (dual.directionOut_<0)
     3275                      effectiveRhs[i]=rowUpper[i];
     3276                    else
     3277                      effectiveRhs[i]=rowLower[i];
     3278                  } else {
     3279                    dual.ray_[i]=0.0;
     3280                  }
     3281                }
     3282              }
     3283            }
     3284            dual.times(-1.0,farkas2,effectiveRhs2);
     3285            double bSum2=0.0;
     3286            for (int i=0;i<numberRows;i++) {
     3287              bSum2 += effectiveRhs2[i]*dual.ray_[i];
     3288            }
     3289            printf("Alternate bSum %g\n",bSum2);
     3290            memset(farkas2,0,numberColumns*sizeof(double));
     3291            dual.transposeTimes(-1.0,dual.ray_,farkas2);
     3292            int nBad=0;
     3293            double largest=-1.0;
     3294            double smallest=1.0e30;
     3295            for (int i=0;i<numberColumns;i++) {
     3296              //if (fabs(farkas[i])>1.0e-12)
     3297              //printf("%d ptr farkas %g dual farkas %g\n",i,farkas[i],farkas2[i]);
     3298              if (fabs(farkas[i]-farkas2[i])>1.0e-7) {
     3299                nBad++;
     3300                largest=CoinMax(largest,fabs(farkas[i]-farkas2[i]));
     3301                smallest=CoinMin(smallest,fabs(farkas[i]-farkas2[i]));
     3302                //printf("%d ptr farkas %g dual farkas %g\n",i,farkas[i],farkas2[i]);
     3303              }
     3304            }
     3305            if (nBad)
     3306              printf("%d farkas difference %g to %g\n",nBad,smallest,largest);
     3307            dual.primal();
     3308            assert (dual.problemStatus_==1);
     3309            assert (!nBad);
     3310          }
     3311        }
     3312        const char * integerInformation = modelPtr_->integerType_;
     3313        assert (integerInformation);
     3314        int * conflict = new int[numberColumns];
     3315        double * sort = new double [numberColumns];
     3316        double relax=0.0;
     3317        int nConflict=0;
     3318        int nOriginal=0;
     3319        int nFixed=0;
     3320        for (int iColumn=0;iColumn<numberColumns;iColumn++) {
     3321          double thisRelax = 0.0;
     3322          if (integerInformation[iColumn]) {
     3323            if ((debugMode&1)!=0)
     3324            printf("%d status %d %g <= %g <=%g (orig %g, %g) farkas %g\n",
     3325                   iColumn,modelPtr_->getStatus(iColumn),columnLower[iColumn],
     3326                   modelPtr_->columnActivity_[iColumn],columnUpper[iColumn],
     3327                   originalLower[iColumn],originalUpper[iColumn],
     3328                   farkas[iColumn]);
     3329            double gap = originalUpper[iColumn]-originalLower[iColumn];
     3330            if (!gap)
     3331              continue;
     3332            if (gap==columnUpper[iColumn]-columnLower[iColumn])
     3333              nOriginal++;
     3334            if (columnUpper[iColumn]==columnLower[iColumn])
     3335              nFixed++;
     3336            if (fabs(farkas[iColumn])<1.0e-15) {
     3337              farkas[iColumn]=0.0;
     3338              continue;
     3339            }
     3340            // temp
     3341            if (gap>=20000.0&&false) {
     3342              // can't use
     3343              if (farkas[iColumn]<0.0) {
     3344                assert(originalLower[iColumn]-columnLower[iColumn]<=0.0);
     3345                // farkas is negative - relax lower bound all way
     3346                relax +=
     3347                  farkas[iColumn]*(originalLower[iColumn]-columnLower[iColumn]);
     3348              } else {
     3349                assert(originalUpper[iColumn]-columnUpper[iColumn]>=0.0);
     3350                // farkas is positive - relax upper bound all way
     3351                relax +=
     3352                  farkas[iColumn]*(originalUpper[iColumn]-columnUpper[iColumn]);
     3353              }
     3354              continue;
     3355            }
     3356            if (originalLower[iColumn]==columnLower[iColumn]) {
     3357              // may need to be careful if odd basic (due to local cut)
     3358              if (farkas[iColumn]>0.0/*&&(modelPtr_->getStatus(iColumn)==ClpSimplex::atUpperBound
     3359                                        ||modelPtr_->getStatus(iColumn)==ClpSimplex::isFixed
     3360                                        ||iColumn==sequenceOut)*/) {
     3361                // farkas is positive - add to list
     3362                gap=originalUpper[iColumn]-columnUpper[iColumn];
     3363                if (gap) {
     3364                  sort[nConflict]=-farkas[iColumn]*gap;
     3365                  conflict[nConflict++]=iColumn;
     3366                }
     3367                //assert (gap>columnUpper[iColumn]-columnLower[iColumn]);
     3368              }
     3369            } else if (originalUpper[iColumn]==columnUpper[iColumn]) {
     3370              // may need to be careful if odd basic (due to local cut)
     3371              if (farkas[iColumn]<0.0/*&&(modelPtr_->getStatus(iColumn)==ClpSimplex::atLowerBound
     3372                                        ||modelPtr_->getStatus(iColumn)==ClpSimplex::isFixed
     3373                                        ||iColumn==sequenceOut)*/) {
     3374                // farkas is negative - add to list
     3375                gap=columnLower[iColumn]-originalLower[iColumn];
     3376                if (gap) {
     3377                  sort[nConflict]=farkas[iColumn]*gap;
     3378                  conflict[nConflict++]=iColumn;
     3379                }
     3380                //assert (gap>columnUpper[iColumn]-columnLower[iColumn]);
     3381              }
     3382            } else {
     3383              // can't use
     3384              if (farkas[iColumn]<0.0) {
     3385                assert(originalLower[iColumn]-columnLower[iColumn]<=0.0);
     3386                // farkas is negative - relax lower bound all way
     3387                thisRelax =
     3388                  farkas[iColumn]*(originalLower[iColumn]-columnLower[iColumn]);
     3389              } else {
     3390                assert(originalUpper[iColumn]-columnUpper[iColumn]>=0.0);
     3391                // farkas is positive - relax upper bound all way
     3392                thisRelax =
     3393                  farkas[iColumn]*(originalUpper[iColumn]-columnUpper[iColumn]);
     3394              }
     3395            }
     3396          } else {
     3397            // not integer - but may have been got at
     3398            double gap = originalUpper[iColumn]-originalLower[iColumn];
     3399            if (gap>columnUpper[iColumn]-columnLower[iColumn]) {
     3400              // can't use
     3401              if (farkas[iColumn]<0.0) {
     3402                assert(originalLower[iColumn]-columnLower[iColumn]<=0.0);
     3403                // farkas is negative - relax lower bound all way
     3404                thisRelax =
     3405                  farkas[iColumn]*(originalLower[iColumn]-columnLower[iColumn]);
     3406              } else {
     3407                assert(originalUpper[iColumn]-columnUpper[iColumn]>=0.0);
     3408                // farkas is positive - relax upper bound all way
     3409                thisRelax =
     3410                  farkas[iColumn]*(originalUpper[iColumn]-columnUpper[iColumn]);
     3411              }
     3412            }
     3413          }
     3414          assert (thisRelax>=0.0);
     3415          relax += thisRelax;
     3416        }
     3417        if (relax+bSum>-1.0e-4||!nConflict) {
     3418          if (relax+bSum>-1.0e-4) {
     3419#if PRINT_CONFLICT>1 //ndef NDEBUG
     3420            printf("General integers relax bSum to %g\n",relax+bSum);
     3421#endif
     3422          } else {
     3423#if PRINT_CONFLICT
     3424            printf("All variables relaxed and still infeasible - what does this mean?\n");
     3425#endif
     3426            int nR=0;
     3427            for (int i=0;i<numberRows;i++) {
     3428              if (fabs(ray[i])>1.0e-10)
     3429                nR++;
     3430              else
     3431                ray[i]=0.0;
     3432            }
     3433            int nC=0;
     3434            for (int i=0;i<numberColumns;i++) {
     3435              if (fabs(farkas[i])>1.0e-10)
     3436                nC++;
     3437              else
     3438                farkas[i]=0.0;
     3439            }
     3440            if (nR<3&&nC<5) {
     3441              printf("BAD %d nonzero rows, %d nonzero columns\n",nR,nC);
     3442            }
     3443          }
     3444        } else if (nConflict < 1000) {
     3445#if PRINT_CONFLICT>1 //ndef NDEBUG
     3446          if (nConflict<5)
     3447            printf("BOUNDS violation bSum %g (relaxed %g) - %d at original bounds, %d fixed - %d in conflict\n",bSum,
     3448                   relax+bSum,nOriginal,nFixed,nConflict);
     3449#endif
     3450          CoinSort_2(sort,sort+nConflict,conflict);
     3451          if ((debugMode&4)!=0) {
     3452            double * temp = new double[numberColumns];
     3453            int * tempC = new int[numberColumns];
     3454            int n=0;
     3455            for (int j=0;j<nConflict;j++) {
     3456              int i=conflict[j];
     3457              if (fabs(farkas[i])>1.0e-12) {
     3458                temp[n]=farkas[i];
     3459                tempC[n++]=i;
     3460              }
     3461            }
     3462            debugModel.addRow(n,tempC,temp,bSum,bSum);
     3463            delete [] tempC;
     3464            delete [] temp;
     3465          }
     3466          int nC=nConflict;
     3467          for (int i=0;i<nConflict;i++) {
     3468            int iColumn=conflict[i];
     3469            if (fabs(sort[i])!=fabs(farkas[iColumn])&&originalUpper[iColumn]==1.0)
     3470              printf("odd %d %g %d %g\n",i,sort[i],iColumn,farkas[iColumn]);
     3471          }
     3472          bSum+=relax;
     3473          double saveBsum = bSum;
     3474          // we can't use same values
     3475          double totalChange=0;
     3476          while (nConflict) {
     3477            double last = -sort[nConflict-1];
     3478            int kConflict=nConflict;
     3479            while (kConflict) {
     3480              double change = -sort[kConflict-1];
     3481              if (change>last+1.0e-5)
     3482                break;
     3483              totalChange += change;
     3484              kConflict--;
     3485            }
     3486            if (bSum+totalChange>-1.0e-4)
     3487              break;
     3488#if 0
     3489            //int iColumn=conflict[nConflict-1];
     3490            double change=-sort[nConflict-1];
     3491            if (bSum+change>-1.0e-4)
     3492              break;
     3493            nConflict--;
     3494            bSum += change;
     3495#else
     3496            nConflict=kConflict;
     3497            bSum += totalChange;
     3498#endif
     3499          }
     3500          if (!nConflict) {
     3501            int nR=0;
     3502            for (int i=0;i<numberRows;i++) {
     3503              if (fabs(ray[i])>1.0e-10)
     3504                nR++;
     3505              else
     3506                ray[i]=0.0;
     3507            }
     3508            int nC=0;
     3509            for (int i=0;i<numberColumns;i++) {
     3510              if (fabs(farkas[i])>1.0e-10)
     3511                nC++;
     3512              else
     3513                farkas[i]=0.0;
     3514            }
     3515#if 1 //PRINT_CONFLICT>1 //ndef NDEBUG
     3516            {
     3517              printf("BAD2 - zero nConflict %d nonzero rows, %d nonzero columns\n",nR,nC);
     3518            }
     3519#endif
     3520          }
     3521          // no point doing if no reduction (or big?) ?
     3522          if (nConflict<nC+1&&nConflict<100&&nConflict) {
     3523            cut=new OsiRowCut();
     3524            cut->setUb(COIN_DBL_MAX);
     3525            if (!typeCut) {
     3526              double lo=1.0;
     3527              for (int i=0;i<nConflict;i++) {
     3528                int iColumn = conflict[i];
     3529                if (originalLower[iColumn]==columnLower[iColumn]) {
     3530                  // must be at least one higher
     3531                  sort[i]=1.0;
     3532                  lo += originalLower[iColumn];
     3533                } else {
     3534                  // must be at least one lower
     3535                  sort[i]=-1.0;
     3536                  lo -= originalUpper[iColumn];
     3537                }
     3538              }
     3539              cut->setLb(lo);
     3540              cut->setRow(nConflict,conflict,sort);
     3541#if PRINT_CONFLICT
     3542              printf("CUT has %d (started at %d) - final bSum %g\n",nConflict,nC,bSum);
     3543#endif
     3544              if ((debugMode&4)!=0) {
     3545                debugModel.addRow(nConflict,conflict,sort,lo,COIN_DBL_MAX);
     3546                debugModel.writeMps("bad.mps");
     3547              }
     3548            } else {
     3549              // just save for use later
     3550              // first take off small
     3551              int nC2=nC;
     3552              while (nC2) {
     3553                //int iColumn=conflict[nConflict-1];
     3554                double change=-sort[nC2-1];
     3555                if (saveBsum+change>-1.0e-4||change>1.0e-4)
     3556                  break;
     3557                nC2--;
     3558                saveBsum += change;
     3559              }
     3560              cut->setLb(saveBsum);
     3561              for (int i=0;i<nC2;i++) {
     3562                int iColumn = conflict[i];
     3563                sort[i]=farkas[iColumn];
     3564              }
     3565              cut->setRow(nC2,conflict,sort);
     3566              // mark as globally valid
     3567              cut->setGloballyValid();
     3568#if PRINT_CONFLICT>1 //ndef NDEBUG
     3569              printf("Stem CUT has %d (greedy %d - with small %d) - saved bSum %g final greedy bSum %g\n",
     3570                     nC2,nConflict,nC,saveBsum,bSum);
     3571#endif
     3572            }
     3573          }
     3574        }
     3575        delete [] conflict;
     3576        delete [] sort;
     3577      }
     3578      delete [] farkas;
     3579    } else {
     3580#if PRINT_CONFLICT>1 //ndef NDEBUG
     3581      printf("Bad dual ray\n");
     3582#endif
     3583    }
     3584    modelPtr_->deleteRay();
     3585  }
     3586  return cut;
     3587}
     3588#endif
     3589
    27073590//#############################################################################
    27083591// Problem information methods (original data)
     
    38564739  return modelPtr_;
    38574740}
    3858 
    38594741//-------------------------------------------------------------------
    38604742
     
    38954777fakeObjective_(NULL)
    38964778{
    3897   //printf("in default %x\n",this);
     4779  //printf("%p %d null constructor\n",this,xxxxxx);xxxxxx++;
    38984780  modelPtr_=NULL;
    38994781  notOwned_=false;
     
    39554837fakeMinInSimplex_(rhs.fakeMinInSimplex_)
    39564838{
    3957   //printf("in copy %x - > %x\n",&rhs,this);
     4839  //printf("%p %d copy constructor %p\n",this,xxxxxx,&rhs);xxxxxx++;
    39584840  if ( rhs.modelPtr_  )
    39594841    modelPtr_ = new ClpSimplex(*rhs.modelPtr_);
     
    40684950OsiClpSolverInterface::~OsiClpSolverInterface ()
    40694951{
    4070   //printf("in destructor %x\n",this);
     4952  //printf("%p destructor\n",this);
    40714953  freeCachedResults();
    40724954  if (!notOwned_)
     
    65447426      }
    65457427      // set normal
    6546       specialOptions_ &= (2047|3*8192|15*65536|2097152|4194304);
     7428      specialOptions_ &= (2047|7*8192|15*65536|2097152|4194304);
    65477429      if (otherInformation!=NULL) {
    65487430        int * array = static_cast<int *> (otherInformation);
     
    67647646                        nBound,moreBounds,tightenBounds);
    67657647#endif
    6766   bool inCbcOrOther = (modelPtr_->specialOptions()&0x03000000)!=0;
     7648  int inCbcOrOther = ((modelPtr_->specialOptions()&0x03000000)!=0) ? 1 : 0;
     7649  if((modelPtr_->specialOptions()&0x01200000)==0x1200000)
     7650    inCbcOrOther=2;
    67677651  if (small) {
    67687652    small->specialOptions_ |= 262144;
     
    68247708    if (problemStatus>=0&&problemStatus<=2) {
    68257709      modelPtr_->setProblemStatus(problemStatus);
    6826       if (!inCbcOrOther||!problemStatus) {
     7710      if (inCbcOrOther!=1||!problemStatus) {
    68277711        // Scaling may have changed - if so pass across
    68287712        if (modelPtr_->scalingFlag()==4)
    68297713          modelPtr_->scaling(small->scalingFlag());
    68307714        static_cast<ClpSimplexOther *> (modelPtr_)->afterCrunch(*small,whichRow,whichColumn,nBound);
    6831         if ((specialOptions_&1048576)==0) {
     7715        if ((specialOptions_&1048576)==0 && !inCbcOrOther) {
    68327716          // get correct rays
    68337717          if (problemStatus==2)
     
    68447728            double * ray = new double [numberRows];
    68457729            memset(ray,0,numberRows*sizeof(double));
     7730            double multiplier=1.0;
     7731            if (small->rowScale_) {
     7732              if (small->sequenceOut_<small->numberColumns_)
     7733                multiplier=small->columnScale_[small->sequenceOut_];
     7734              else
     7735                multiplier=small->rowScale_[small->sequenceOut_-small->numberColumns_];
     7736            }
    68467737            for (int i = 0; i < numberRows2; i++) {
     7738              small->ray_[i] *= multiplier;
    68477739              int iRow = whichRow[i];
    68487740              ray[iRow] = small->ray_[i];
     7741            }
     7742            {
     7743              // Column copy of matrix
     7744              const double * element = small->matrix()->getElements();
     7745              const int * row = small->matrix()->getIndices();
     7746              const CoinBigIndex * columnStart = small->matrix()->getVectorStarts();
     7747              const int * columnLength = small->matrix()->getVectorLengths();
     7748              int n=0,k=0,nn=0;
     7749              for (int i=0;i<small->numberColumns_;i++) {
     7750                double sum = 0.0;
     7751                for (CoinBigIndex j = columnStart[i];
     7752                     j < columnStart[i] + columnLength[i]; j++) {
     7753                  sum += small->ray_[row[j]]*element[j];
     7754                }
     7755                if (fabs(sum)>1.0e-7) {
     7756                  if (small->getColumnStatus(i) == ClpSimplex::basic) {
     7757                    n++;
     7758                    if (fabs(1.0-fabs(sum))>1.0e-7)
     7759                      nn++;
     7760                  } else if (fabs(sum)>1.0e-7) {
     7761                    k++;
     7762                  }
     7763                  //printf("%d %g\n",i,sum);
     7764                }
     7765              }
     7766              //printf("small %d basic (%d non-unit) %d non-basic\n",n,nn,k);
    68497767            }
    68507768            // Column copy of matrix
     
    68547772            const int * columnLength = getMatrixByCol()->getVectorLengths();
    68557773            // translate
    6856             //pivotRow=whichRow[pivotRow];
    6857             //modelPtr_->spareIntArray_[3]=pivotRow;
    6858             int pivotRow=-1;
    68597774            for (int jRow = nBound; jRow < 2 * numberRows; jRow++) {
    68607775              int iRow = whichRow[jRow];
     
    68717786                  }
    68727787                }
    6873                 if (iRow!=pivotRow) {
    6874                   ray[iRow] = -sum / value;
    6875                 } else {
    6876                   printf("what now - direction %d wanted %g sum %g value %g\n",
    6877                          small->directionOut_,ray[iRow],
    6878                          sum,value);
    6879                 }
     7788                ray[iRow] = -sum / value;
    68807789              }
    68817790            }
     7791            int n=0,k=0,nn=0;
    68827792            for (int i=0;i<modelPtr_->numberColumns_;i++) {
     7793              double sum = 0.0;
     7794              for (CoinBigIndex j = columnStart[i];
     7795                   j < columnStart[i] + columnLength[i]; j++) {
     7796                sum += ray[row[j]]*element[j];
     7797              }
     7798              if (modelPtr_->getStatus(i)==ClpSimplex::basic&&
     7799                  fabs(sum)>1.0e-7) {
     7800                n++;
     7801                if (fabs(1.0-fabs(sum))>1.0e-7)
     7802                  nn++;
     7803              } else if (fabs(sum)>1.0e-7) {
     7804                k++;
     7805              }
    68837806              if (modelPtr_->getStatus(i)!=ClpSimplex::basic&&
    68847807                  modelPtr_->columnLower_[i]==modelPtr_->columnUpper_[i])
    68857808                modelPtr_->setStatus(i,ClpSimplex::isFixed);
    68867809            }
     7810            //printf("%d basic (%d non-unit) %d non-basic\n",n,nn,k);
    68877811            modelPtr_->ray_=ray;
    68887812            modelPtr_->directionOut_=small->directionOut_;
     7813            if (small->sequenceOut_<small->numberColumns_)
     7814              modelPtr_->sequenceOut_ = whichColumn[small->sequenceOut_];
     7815            else
     7816              modelPtr_->sequenceOut_ = whichRow[small->sequenceOut_-small->numberColumns_]
     7817                + modelPtr_->numberColumns_;
    68897818          }
    68907819        }
     
    79238852            }
    79248853            // See if can do anything
    7925             {
     8854            if (1) {
    79268855              /*
    79278856                If kNode is on second branch then
     
    79428871                  assert (node.way_==2||node.way_==-2);
    79438872                  // we don't know which branch it was - look at current bounds
    7944                   if (upper[iColumn]<value&&node.lower_[variable]<upper[iColumn]) {
     8873                  if (upper[iColumn]<value&&node.upper_[variable]>upper[iColumn]) {
    79458874                    // must have been down branch
    79468875                    if (djValue>1.0e-3||solution[iColumn]<upper[iColumn]-1.0e-5) {
     
    80038932                    }
    80048933                    // we don't know which branch it was - look at current bounds
    8005                   } else if (lower[iColumn]>value&&node.upper_[variable]>lower[iColumn]) {
     8934                  } else if (lower[iColumn]>value&&node.lower_[variable]<lower[iColumn]) {
    80068935                    // must have been up branch
    80078936                    if (djValue<-1.0e-3||solution[iColumn]>lower[iColumn]+1.0e-5) {
     
    80749003            //assert(!getIterationCount());
    80759004            if ((numberNodes%1000)==0)
    8076               printf("%d nodes, redundant down %d (%d) up %d (%d) tree size %d\n",
    8077                      numberNodes,nRedundantDown,nRedundantDown2,nRedundantUp,nRedundantUp2,branchingTree.size());
     9005              printf("%d nodes, redundant down %d (%d) up %d (%d) tree size %d best solution %g\n",
     9006                     numberNodes,nRedundantDown,nRedundantDown2,nRedundantUp,nRedundantUp2,
     9007                     branchingTree.size(),bestNode.objectiveValue_);
    80789008#else
    80799009            if ((numberNodes%1000)==0)
     
    81019031          } else {
    81029032            // infeasible
     9033            // could use conflict cut code to find redundant
    81039034            delete ws;
    81049035          }
     
    81579088OsiClpSolverInterface::setSpecialOptions(unsigned int value)
    81589089{
     9090#ifdef NO_CRUNCH2
     9091  value &= ~1;
     9092#endif
    81599093  if ((value&131072)!=0&&(specialOptions_&131072)==0) {
    81609094    // Try and keep scaling factors around
     
    82759209OsiClpSolverInterface::tightenBounds(int lightweight)
    82769210{
     9211#if FUNNY_BRANCHING3>100
     9212  return 0;
     9213#endif
    82779214  if (!integerInformation_||(specialOptions_&262144)!=0)
    82789215    return 0; // no integers
     
    88729809  return nTightened;
    88739810}
     9811// See if any integer variables make infeasible other way
     9812int
     9813OsiClpSolverInterface::infeasibleOtherWay(char * whichWay)
     9814{
     9815  //CoinPackedMatrix matrixByRow(*getMatrixByRow());
     9816  int numberRows = getNumRows();
     9817  int numberColumns = getNumCols();
     9818
     9819  int iRow,iColumn;
     9820
     9821  // Row copy
     9822  //const double * elementByRow = matrixByRow.getElements();
     9823  //const int * column = matrixByRow.getIndices();
     9824  //const CoinBigIndex * rowStart = matrixByRow.getVectorStarts();
     9825  //const int * rowLength = matrixByRow.getVectorLengths();
     9826
     9827  const double * columnUpper = getColUpper();
     9828  const double * columnLower = getColLower();
     9829  const double * rowUpper = getRowUpper();
     9830  const double * rowLower = getRowLower();
     9831
     9832  // Column copy of matrix
     9833  const double * element = getMatrixByCol()->getElements();
     9834  const int * row = getMatrixByCol()->getIndices();
     9835  const CoinBigIndex * columnStart = getMatrixByCol()->getVectorStarts();
     9836  const int * columnLength = getMatrixByCol()->getVectorLengths();
     9837  const double *objective = getObjCoefficients() ;
     9838  double direction = getObjSense();
     9839  double tolerance = 1.0e-6;
     9840  double * down = new double [3*numberRows];
     9841  double * up = down + numberRows;
     9842  double * sum = up + numberRows;
     9843  int * nLowerInf = new int [2*numberRows];
     9844  int * nUpperInf = nLowerInf + numberRows;
     9845  CoinZeroN(down,numberRows);
     9846  CoinZeroN(up,numberRows);
     9847  CoinZeroN(sum,numberRows);
     9848  for (int i=0;i<numberRows;i++) {
     9849    nLowerInf[i]=-1;
     9850    nUpperInf[i]=-1;
     9851  }
     9852  double infinity = getInfinity();
     9853  for (iColumn=0;iColumn<numberColumns;iColumn++) {
     9854    CoinBigIndex start = columnStart[iColumn];
     9855    CoinBigIndex end = start + columnLength[iColumn];
     9856    double lower = columnLower[iColumn];
     9857    double upper = columnUpper[iColumn];
     9858    if (lower==upper) {
     9859      for (CoinBigIndex j=start;j<end;j++) {
     9860        int iRow = row[j];
     9861        double value = element[j];
     9862        down[iRow] += value*lower;
     9863        up[iRow] += value*lower;
     9864        sum[iRow] += 2.0*fabs(value*lower);
     9865      }
     9866    } else {
     9867      for (CoinBigIndex j=start;j<end;j++) {
     9868        int iRow = row[j];
     9869        double value = element[j];
     9870        if (value>0.0) {
     9871          if (lower!=-infinity) {
     9872            down[iRow] += value*lower;
     9873            sum[iRow]+=fabs(value*lower);
     9874          } else if (nLowerInf[iRow]==-1) {
     9875            nLowerInf[iRow]=iColumn;
     9876          } else {
     9877            nLowerInf[iRow]=-2;
     9878          }
     9879          if (upper!=infinity) {
     9880            up[iRow] += value*upper;
     9881            sum[iRow]+=fabs(value*upper);
     9882          } else if (nUpperInf[iRow]==-1) {
     9883            nUpperInf[iRow]=iColumn;
     9884          } else {
     9885            nUpperInf[iRow]=-2;
     9886          }
     9887        } else {
     9888          if (upper!=infinity) {
     9889            down[iRow] += value*upper;
     9890            sum[iRow]+=fabs(value*upper);
     9891          } else if (nLowerInf[iRow]==-1) {
     9892            nLowerInf[iRow]=iColumn;
     9893          } else {
     9894            nLowerInf[iRow]=-2;
     9895          }
     9896          if (lower!=-infinity) {
     9897            up[iRow] += value*lower;
     9898            sum[iRow]+=fabs(value*lower);
     9899          } else if (nUpperInf[iRow]==-1) {
     9900            nUpperInf[iRow]=iColumn;
     9901          } else {
     9902            nUpperInf[iRow]=-2;
     9903          }
     9904        }
     9905      }
     9906    }
     9907  }
     9908  int nOtherWay=0;
     9909  for (iColumn=0;iColumn<numberColumns;iColumn++) {
     9910    char upDown = whichWay[iColumn];
     9911    if (!upDown)
     9912      continue;
     9913    // take off contribution of column
     9914    double lower = columnLower[iColumn];
     9915    double upper = columnUpper[iColumn];
     9916    double newLower;
     9917    double newUpper;
     9918    if (upDown==1) {
     9919      newLower=-infinity;
     9920      newUpper=lower -1.0;
     9921    } else {
     9922      newLower=upper+1.0;
     9923      newUpper=infinity;
     9924    }
     9925    CoinBigIndex start = columnStart[iColumn];
     9926    CoinBigIndex end = start + columnLength[iColumn];
     9927    if (lower<-1.0e8&&upper>1.0e8)
     9928      continue; // Could do severe damage to accuracy
     9929    for (CoinBigIndex j=start;j<end;j++) {
     9930      int iRow = row[j];
     9931      double value = element[j];
     9932      double thisDown = down[iRow];
     9933      double thisUp = up[iRow];
     9934      if (value>0.0) {
     9935        if (nLowerInf[iRow]==-1)
     9936          thisDown -= value*lower;
     9937        else if (nLowerInf[iRow]==-2||
     9938                 (nLowerInf[iRow]!=iColumn))
     9939          thisDown = -infinity;
     9940        if (nUpperInf[iRow]==-1)
     9941          thisUp -= value*upper;
     9942        else if (nUpperInf[iRow]==-2||
     9943                 (nUpperInf[iRow]!=iColumn))
     9944          thisUp = infinity;
     9945        thisDown += newLower*value;
     9946        thisUp += newUpper*value;
     9947      } else {
     9948        if (nLowerInf[iRow]==-1)
     9949          thisUp -= value*lower;
     9950        else if (nLowerInf[iRow]==-2||
     9951                 (nLowerInf[iRow]!=iColumn))
     9952          thisUp = infinity;
     9953        if (nUpperInf[iRow]==-1)
     9954          thisDown -= value*upper;
     9955        else if (nUpperInf[iRow]==-2||
     9956                 (nUpperInf[iRow]!=iColumn))
     9957          thisDown = -infinity;
     9958        thisDown += newUpper*value;
     9959        thisUp += newLower*value;
     9960      }
     9961      if (thisDown>rowUpper[iRow]) {
     9962        if (thisDown>rowUpper[iRow]+tolerance+1.0e-8*sum[iRow]) {
     9963          nOtherWay++;
     9964          whichWay[iColumn] *=2;
     9965        }
     9966      } else if (thisUp<rowLower[iRow]) {
     9967        if (thisUp<rowUpper[iRow]-tolerance-1.0e-8*sum[iRow]) {
     9968          nOtherWay++;
     9969          whichWay[iColumn] *=2;
     9970        }
     9971      }
     9972    }
     9973  }
     9974  delete [] nLowerInf;
     9975  delete [] down;
     9976  return nOtherWay;
     9977}
    88749978// Return number of entries in L part of current factorization
    88759979CoinBigIndex
     
    941010514      return false;
    941110515    } else if (phase_<2) {
    9412       if (model_->numberIterations()> model_->baseIteration()+2*model_->numberRows()+model_->numberColumns()+2000||
     10516      if (model_->numberIterations()> model_->baseIteration()+2*model_->numberRows()+model_->numberColumns()+100000||
    941310517          model_->largestDualError()>=1.0e-1) {
    941410518        // had model_->numberDualInfeasibilitiesWithoutFree()||
     
    944110545    } else {
    944210546      assert (phase_==2);
    9443       if (model_->numberIterations()> model_->baseIteration()+3*model_->numberRows()+model_->numberColumns()+2000||
     10547      if (model_->numberIterations()> model_->baseIteration()+3*model_->numberRows()+model_->numberColumns()+100000||
    944410548          model_->largestPrimalError()>=1.0e3) {
    944510549#ifdef COIN_DEVELOP
     
    945310557  } else {
    945410558    // primal
    9455     if (model_->numberIterations()<model_->baseIteration()+2*model_->numberRows()+model_->numberColumns()+4000) {
     10559    if (model_->numberIterations()<model_->baseIteration()+2*model_->numberRows()+model_->numberColumns()+100000) {
    945610560      return false;
    945710561    } else if (phase_<2) {
    9458       if (model_->numberIterations()> model_->baseIteration()+3*model_->numberRows()+2000+
     10562      if (model_->numberIterations()> model_->baseIteration()+3*model_->numberRows()+20000+
    945910563          model_->numberColumns()&&
    946010564          model_->numberDualInfeasibilitiesWithoutFree()>0&&
     
    947010574    } else {
    947110575      assert (phase_==2);
    9472       if (model_->numberIterations()> model_->baseIteration()+3*model_->numberRows()+2000||
     10576      if (model_->numberIterations()> model_->baseIteration()+3*model_->numberRows()+20000||
    947310577          model_->largestPrimalError()>=1.0e3) {
    947410578#ifdef COIN_DEVELOP
     
    948610590{
    948710591  inTrouble_=true;
     10592}
     10593// are we in trouble
     10594bool
     10595OsiClpDisasterHandler::inTrouble() const
     10596{
     10597  return inTrouble_|(osiModel_->getModelPtr()->problemStatus()==4);
    948810598}
    948910599// Type of disaster 0 can fix, 1 abort
  • trunk/Clp/src/OsiClp/OsiClpSolverInterface.hpp

    r1972 r2070  
    340340  /// Sets integer tolerance and increment
    341341  void setStuff(double tolerance,double increment);
     342  /// Return a conflict analysis cut from small model
     343  OsiRowCut * smallModelCut(const double * originalLower, const double * originalUpper,
     344                            int numberRowsAtContinuous,const int * whichGenerator,
     345                            int typeCut=0);
     346  /** Return a conflict analysis cut from model
     347      If type is 0 then genuine cut, if 1 then only partially processed
     348   */
     349  OsiRowCut * modelCut(const double * originalLower, const double * originalUpper,
     350                       int numberRowsAtContinuous,const int * whichGenerator,
     351                       int typeCut=0);
    342352  //@}
    343353 
     
    10801090  */
    10811091  virtual int tightenBounds(int lightweight=0);
     1092  /// See if any integer variables make infeasible other way
     1093  int infeasibleOtherWay(char * whichWay);
    10821094  /// Return number of entries in L part of current factorization
    10831095  virtual CoinBigIndex getSizeL() const;
     
    14591471  { return phase_;}
    14601472  /// are we in trouble
    1461   inline bool inTrouble() const
    1462   { return inTrouble_;}
     1473  bool inTrouble() const;
    14631474 
    14641475  //@}
Note: See TracChangeset for help on using the changeset viewer.