Changeset 2025 for trunk


Ignore:
Timestamp:
Mar 19, 2014 8:49:55 AM (6 years ago)
Author:
forrest
Message:

mainly dantzig-wolfe - also Clp strong branching (rather than OsiClp?)

Location:
trunk/Clp/src
Files:
13 edited

Legend:

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

    r1910 r2025  
    484484     }
    485485     if (model_->factorization()->pivots() && model_->numberPrimalInfeasibilities())
    486           tolerance = CoinMax(tolerance, 1.0e-10 * model_->infeasibilityCost());
     486          tolerance = CoinMax(tolerance, 1.0e-15 * model_->infeasibilityCost());
    487487     tolerance *= tolerance; // as we are using squares
    488488
     
    20322032     }
    20332033     if (model_->factorization()->pivots() && model_->numberPrimalInfeasibilities())
    2034        tolerance = CoinMax(tolerance, 1.0e-10 * model_->infeasibilityCost());
     2034       tolerance = CoinMax(tolerance, 1.0e-15 * model_->infeasibilityCost());
    20352035     // So partial pricing can use
    20362036     model_->setCurrentDualTolerance(tolerance);
  • trunk/Clp/src/AbcSimplex.cpp

    r2024 r2025  
    28362836  double relaxedTolerance = primalTolerance_;
    28372837  // we can't really trust infeasibilities if there is primal error
    2838   double error = CoinMin(1.0e-2, largestPrimalError_);
     2838  double error = CoinMin(1.0e-2, CoinMax(largestPrimalError_,5.0*primalTolerance_));
    28392839  // allow tolerance at least slightly bigger than standard
    28402840  relaxedTolerance = relaxedTolerance +  error;
     
    28872887  bestPossibleImprovement_ = 0.0;
    28882888  // we can't really trust infeasibilities if there is dual error
    2889   double error = CoinMin(1.0e-2, largestDualError_);
     2889  double error = CoinMin(1.0e-2, CoinMax(largestDualError_,5.0*dualTolerance_));
    28902890  // allow tolerance at least slightly bigger than standard
    28912891  double relaxedTolerance = currentDualTolerance_ +  error;
  • trunk/Clp/src/ClpModel.cpp

    r2024 r2025  
    14271427          if (status_) {
    14281428               if (numberColumns_ + newSize) {
     1429#define CLP_TIDY_DELETE_ROWS
     1430#ifdef CLP_TIDY_DELETE_ROWS
     1431                    // try and get right number of basic
     1432                    int nChange=0;
     1433                    unsigned char * rowStatus = status_+numberColumns_;
     1434                    for (int i=0;i<number;i++) {
     1435                      int iRow = which[i];
     1436                      if ((rowStatus[iRow]&7) !=1)
     1437                        nChange++;
     1438                    }
     1439                    // take out slacks at bound
     1440                    for (int iRow=0;iRow<numberRows_;iRow++) {
     1441                        if (!nChange)
     1442                          break;
     1443                        if ((rowStatus[iRow]&7) ==1) {
     1444                          if (fabs(rowActivity_[iRow]-rowLower_[iRow])<1.0e-8) {
     1445                            rowStatus[iRow]=3;
     1446                            nChange--;
     1447                          } else if (fabs(rowActivity_[iRow]-rowUpper_[iRow])<1.0e-8) {
     1448                            rowStatus[iRow]=2;
     1449                            nChange--;
     1450                          }
     1451                        }
     1452                      }
     1453#endif
    14291454                    unsigned char * tempR  = reinterpret_cast<unsigned char *>
    14301455                                             (deleteChar(reinterpret_cast<char *>(status_) + numberColumns_,
  • trunk/Clp/src/ClpPrimalColumnSteepest.cpp

    r1956 r2025  
    512512     }
    513513     if (model_->factorization()->pivots() && model_->numberPrimalInfeasibilities())
    514           tolerance = CoinMax(tolerance, 1.0e-10 * model_->infeasibilityCost());
     514          tolerance = CoinMax(tolerance, 1.0e-15 * model_->infeasibilityCost());
    515515     tolerance *= tolerance; // as we are using squares
    516516
     
    36983698     }
    36993699     if (model_->factorization()->pivots() && model_->numberPrimalInfeasibilities())
    3700           tolerance = CoinMax(tolerance, 1.0e-10 * model_->infeasibilityCost());
     3700          tolerance = CoinMax(tolerance, 1.0e-15 * model_->infeasibilityCost());
    37013701     // So partial pricing can use
    37023702     model_->setCurrentDualTolerance(tolerance);
  • trunk/Clp/src/ClpSimplex.cpp

    r2024 r2025  
    46164616     dualRowPivot_->setModel(this);
    46174617}
    4618 // Sets row pivot choice algorithm in dual
     4618// Sets column pivot choice algorithm in primal
    46194619void
    46204620ClpSimplex::setPrimalColumnPivotAlgorithm(ClpPrimalColumnPivot & choice)
     
    62646264     model2->setPerturbation(savePerturbation);
    62656265     return model2->status();
     6266}
     6267typedef struct {
     6268     char * spareArrays_;
     6269     ClpFactorization * factorization_;
     6270     int logLevel_;
     6271} ClpHotSaveData;
     6272// Create a hotstart point of the optimization process
     6273void
     6274ClpSimplex::markHotStart(void * &saveStuff)
     6275{
     6276  ClpHotSaveData * saveData = new ClpHotSaveData;
     6277  saveStuff = saveData;
     6278  setProblemStatus(0);
     6279  saveData->logLevel_=logLevel();
     6280  if (logLevel()<2)
     6281    setLogLevel(0);
     6282  // Get space for strong branching
     6283  int size = static_cast<int>((1+4*(numberRows_+numberColumns_))*sizeof(double));
     6284  // and for save of original column bounds
     6285  size += static_cast<int>(2*numberColumns_*sizeof(double));
     6286  size += static_cast<int>((1+4*numberRows_+2*numberColumns_)*sizeof(int));
     6287  size += numberRows_+numberColumns_;
     6288  saveData->spareArrays_ = new char[size];
     6289  // Setup for strong branching
     6290  saveData->factorization_ =
     6291    static_cast<ClpSimplexDual *>(this)->setupForStrongBranching(saveData->spareArrays_,numberRows_,numberColumns_,true);
     6292  double * arrayD = reinterpret_cast<double *> (saveData->spareArrays_);
     6293  arrayD[0]=objectiveValue()* optimizationDirection();
     6294  double * saveSolution = arrayD+1;
     6295  double * saveLower = saveSolution + (numberRows_+numberColumns_);
     6296  double * saveUpper = saveLower + (numberRows_+numberColumns_);
     6297  double * saveObjective = saveUpper + (numberRows_+numberColumns_);
     6298  double * saveLowerOriginal = saveObjective + (numberRows_+numberColumns_);
     6299  double * saveUpperOriginal = saveLowerOriginal + numberColumns_;
     6300  CoinMemcpyN( columnLower(),numberColumns_, saveLowerOriginal);
     6301  CoinMemcpyN( columnUpper(),numberColumns_, saveUpperOriginal);
     6302}
     6303// Optimize starting from the hotstart
     6304void
     6305ClpSimplex::solveFromHotStart(void * saveStuff)
     6306{
     6307  ClpHotSaveData * saveData = reinterpret_cast<ClpHotSaveData *>(saveStuff);
     6308  int iterationLimit = intParam_[ClpMaxNumIteration];
     6309  intParam_[ClpMaxNumIteration] = intParam_[ClpMaxNumIterationHotStart];
     6310  double * arrayD = reinterpret_cast<double *> (saveData->spareArrays_);
     6311  double saveObjectiveValue = arrayD[0];
     6312  double * saveSolution = arrayD+1;
     6313  int number = numberRows_+numberColumns_;
     6314  CoinMemcpyN(saveSolution,number,solutionRegion());
     6315  double * saveLower = saveSolution + (numberRows_+numberColumns_);
     6316  CoinMemcpyN(saveLower,number,lowerRegion());
     6317  double * saveUpper = saveLower + (numberRows_+numberColumns_);
     6318  CoinMemcpyN(saveUpper,number,upperRegion());
     6319  double * saveObjective = saveUpper + (numberRows_+numberColumns_);
     6320  CoinMemcpyN(saveObjective,number,costRegion());
     6321  double * saveLowerOriginal = saveObjective + (numberRows_+numberColumns_);
     6322  double * saveUpperOriginal = saveLowerOriginal + numberColumns_;
     6323  arrayD = saveUpperOriginal + numberColumns_;
     6324  int * savePivot = reinterpret_cast<int *> (arrayD);
     6325  CoinMemcpyN(savePivot,numberRows_,pivotVariable());
     6326  int * whichRow = savePivot+numberRows_;
     6327  int * whichColumn = whichRow + 3*numberRows_;
     6328  int * arrayI = whichColumn + 2*numberColumns_;
     6329  unsigned char * saveStatus = reinterpret_cast<unsigned char *> (arrayI+1);
     6330  CoinMemcpyN(saveStatus,number,statusArray());
     6331  setFactorization(*saveData->factorization_);
     6332  // make sure whatsChanged_ has 1 set
     6333  setWhatsChanged(511);
     6334  double * lowerInternal = lowerRegion();
     6335  double * upperInternal = upperRegion();
     6336  double rhsScale = this->rhsScale();
     6337  const double * columnScale = NULL;
     6338  // and do bounds in case dual needs them
     6339  int iColumn;
     6340  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
     6341    if (columnLower_[iColumn]>saveLowerOriginal[iColumn]) {
     6342      double value = columnLower_[iColumn];
     6343      value *= rhsScale;
     6344      if (columnScale_)
     6345        value /= columnScale_[iColumn];
     6346      lowerInternal[iColumn]=value;
     6347    }
     6348    if (columnUpper_[iColumn]<saveUpperOriginal[iColumn]) {
     6349      double value = columnUpper_[iColumn];
     6350      value *= rhsScale;
     6351      if (columnScale_)
     6352        value /= columnScale_[iColumn];
     6353      upperInternal[iColumn]=value;
     6354    }
     6355  }
     6356  //#define REDO_STEEPEST_EDGE
     6357#ifdef REDO_STEEPEST_EDGE
     6358    if (dynamic_cast<ClpDualRowSteepest *>(dualRowPivot_)) {
     6359      // new copy (should think about keeping around)
     6360      ClpDualRowSteepest steepest;
     6361      setDualRowPivotAlgorithm(steepest);
     6362    }     
     6363#endif
     6364  // Start of fast iterations
     6365  bool alwaysFinish = true; //((specialOptions_&32)==0) ? true : false;
     6366  int saveNumberFake = (static_cast<ClpSimplexDual *>(this))->numberFake_;
     6367  int status = (static_cast<ClpSimplexDual *>(this))->fastDual(alwaysFinish);
     6368  (static_cast<ClpSimplexDual *>(this))->numberFake_ = saveNumberFake;
     6369 
     6370  int probStatus = problemStatus();
     6371  double objValue =objectiveValue() * optimizationDirection();
     6372  CoinAssert (probStatus||objValue<1.0e50);
     6373  // make sure plausible
     6374  double obj = CoinMax(objValue,saveObjectiveValue);
     6375  if (status==10||status<0) {
     6376    // was trying to clean up or something odd
     6377    status=1;
     6378  }
     6379  if (status) {
     6380    // not finished - might be optimal
     6381    checkPrimalSolution(solutionRegion(0),
     6382                              solutionRegion(1));
     6383    objValue =objectiveValue() *
     6384      optimizationDirection();
     6385    obj = CoinMax(objValue,saveObjectiveValue);
     6386    if (!numberDualInfeasibilities()) {
     6387      double limit = 0.0;
     6388      getDblParam(ClpDualObjectiveLimit, limit);
     6389      if (secondaryStatus()==1&&!probStatus&&obj<limit) {
     6390        obj=limit;
     6391        probStatus=3;
     6392      }
     6393      if (!numberPrimalInfeasibilities()&&obj<limit) {
     6394        probStatus=0;
     6395      } else if (probStatus==10) {
     6396        probStatus=3;
     6397      } else if (!numberPrimalInfeasibilities()) {
     6398        probStatus=1; // infeasible
     6399      }
     6400    } else {
     6401      // can't say much
     6402      probStatus=3;
     6403    }
     6404  } else if (!probStatus) {
     6405    if (isDualObjectiveLimitReached())
     6406      probStatus=1; // infeasible
     6407  }
     6408  if (status&&!probStatus) {
     6409    probStatus=3; // can't be sure
     6410  }
     6411  if (probStatus<0)
     6412    probStatus=3;
     6413  setProblemStatus(probStatus);
     6414  setObjectiveValue(obj*optimizationDirection());
     6415  double * solution = primalColumnSolution();
     6416  const double * solution2 = solutionRegion();
     6417  // could just do changed bounds - also try double size scale so can use * not /
     6418  if (!columnScale) {
     6419    for (iColumn=0;iColumn<numberColumns_;iColumn++) {
     6420      solution[iColumn]= solution2[iColumn];
     6421    }
     6422  } else {
     6423    for (iColumn=0;iColumn<numberColumns_;iColumn++) {
     6424      solution[iColumn]= solution2[iColumn]*columnScale[iColumn];
     6425    }
     6426  }
     6427  CoinMemcpyN(saveLowerOriginal,numberColumns_,columnLower_);
     6428  CoinMemcpyN(saveUpperOriginal,numberColumns_,columnUpper_);
     6429  // and back bounds
     6430  CoinMemcpyN(saveLower,number,lowerRegion());
     6431  CoinMemcpyN(saveUpper,number,upperRegion());
     6432  intParam_[ClpMaxNumIteration] = iterationLimit;
     6433}
     6434// Delete the snapshot
     6435void
     6436ClpSimplex::unmarkHotStart(void * saveStuff)
     6437{
     6438  ClpHotSaveData * saveData = reinterpret_cast<ClpHotSaveData *>(saveStuff);
     6439  setLogLevel(saveData->logLevel_);
     6440  deleteRim(0);
     6441  delete saveData->factorization_;
     6442  delete [] saveData->spareArrays_;
     6443  delete saveData;
    62666444}
    62676445/* For strong branching.  On input lower and upper are new bounds
  • trunk/Clp/src/ClpSimplex.hpp

    r2015 r2025  
    306306      startup 1 for values pass
    307307      interrupt whether to pass across interrupt handler
     308      add 10 to return AbcSimplex
    308309  */
    309   void dealWithAbc(int solveType,int startUp,bool interrupt=false);
     310  AbcSimplex * dealWithAbc(int solveType,int startUp,bool interrupt=false);
     311  //void dealWithAbc(int solveType,int startUp,bool interrupt=false);
    310312#endif
    311313     /** This loads a model from a CoinStructuredModel object - returns number of errors.
     
    421423     /// After mini presolve
    422424     void miniPostsolve(const ClpSimplex * presolvedModel,void * info);
    423      /// Scale real objective
    424      void scaleRealObjective(double multiplier);
    425      /// Scale so no RHS (abs not infinite) > value
    426      void scaleRealRhs(double maxValue, double killIfSmaller);
     425     /// mini presolve and solve
     426     void miniSolve(char * rowType, char *columnType,int algorithm, int startUp);
    427427     /** Write the basis in MPS format to the specified file.
    428428         If writeValues true writes values of structurals
     
    484484     /// Sets column pivot choice algorithm in primal
    485485     void setPrimalColumnPivotAlgorithm(ClpPrimalColumnPivot & choice);
     486     /// Create a hotstart point of the optimization process
     487     void markHotStart(void * & saveStuff);
     488     /// Optimize starting from the hotstart
     489     void solveFromHotStart(void * saveStuff);
     490     /// Delete the snapshot
     491     void unmarkHotStart(void * saveStuff);
    486492     /** For strong branching.  On input lower and upper are new bounds
    487493         while on output they are change in objective function values
     
    821827     double scaleObjective(double value);
    822828     /// Solve using Dantzig-Wolfe decomposition and maybe in parallel
    823      int solveDW(CoinStructuredModel * model);
     829    int solveDW(CoinStructuredModel * model, ClpSolve & options);
    824830     /// Solve using Benders decomposition and maybe in parallel
    825      int solveBenders(CoinStructuredModel * model);
     831     int solveBenders(CoinStructuredModel * model, ClpSolve & options);
    826832public:
    827833     /** For advanced use.  When doing iterative solves things can get
     
    12731279     inline bool active(int iRow) const {
    12741280          return ((status_[iRow] & 128) != 0);
     1281     }
     1282     /// To say perturbed
     1283     inline void setPerturbed( int iSequence) {
     1284          status_[iSequence] = static_cast<unsigned char>(status_[iSequence] | 128);
     1285     }
     1286     inline void clearPerturbed( int iSequence) {
     1287          status_[iSequence] = static_cast<unsigned char>(status_[iSequence] & ~128);
     1288     }
     1289     inline bool perturbed(int iSequence) const {
     1290          return ((status_[iSequence] & 128) != 0);
    12751291     }
    12761292     /** Set up status array (can be used by OsiClp).
     
    16651681#define CLP_ABC_BEEN_FEASIBLE 65536
    16661682  int abcState_;
     1683  /// Number of degenerate pivots since last perturbed
     1684  int numberDegeneratePivots_;
    16671685public:
    16681686     /// Spare int array for passing information [0]!=0 switches on
     
    16731691     /// Allow OsiClp certain perks
    16741692     friend class OsiClpSolverInterface;
     1693     /// And OsiCLP
     1694     friend class OsiCLPSolverInterface;
    16751695     //@}
    16761696};
  • trunk/Clp/src/ClpSimplexDual.cpp

    r2015 r2025  
    60456045               status = 1;
    60466046          }
    6047 
    6048           if (scalingFlag_ <= 0) {
    6049                CoinMemcpyN(solution_, numberColumns_, outputSolution[iSolution]);
    6050           } else {
    6051                int j;
    6052                double * sol = outputSolution[iSolution];
    6053                for (j = 0; j < numberColumns_; j++)
    6054                     sol[j] = solution_[j] * columnScale_[j];
    6055           }
     6047          if (outputSolution) {
     6048            if (scalingFlag_ <= 0) {
     6049              CoinMemcpyN(solution_, numberColumns_, outputSolution[iSolution]);
     6050            } else {
     6051              int j;
     6052              double * sol = outputSolution[iSolution];
     6053              for (j = 0; j < numberColumns_; j++)
     6054                sol[j] = solution_[j] * columnScale_[j];
     6055            }
     6056          }
    60566057          outputStatus[iSolution] = status;
    60576058          outputIterations[iSolution] = numberIterations_;
     
    61166117               status = 1;
    61176118          }
    6118           if (scalingFlag_ <= 0) {
    6119                CoinMemcpyN(solution_, numberColumns_, outputSolution[iSolution]);
    6120           } else {
    6121                int j;
    6122                double * sol = outputSolution[iSolution];
    6123                for (j = 0; j < numberColumns_; j++)
    6124                     sol[j] = solution_[j] * columnScale_[j];
    6125           }
     6119          if (outputSolution) {
     6120            if (scalingFlag_ <= 0) {
     6121              CoinMemcpyN(solution_, numberColumns_, outputSolution[iSolution]);
     6122            } else {
     6123              int j;
     6124              double * sol = outputSolution[iSolution];
     6125              for (j = 0; j < numberColumns_; j++)
     6126                sol[j] = solution_[j] * columnScale_[j];
     6127            }
     6128          }
    61266129          outputStatus[iSolution] = status;
    61276130          outputIterations[iSolution] = numberIterations_;
  • trunk/Clp/src/ClpSimplexNonlinear.cpp

    r1929 r2025  
    1010#include "ClpHelperFunctions.hpp"
    1111#include "ClpSimplexNonlinear.hpp"
     12#include "ClpSimplexOther.hpp"
     13#include "ClpSimplexDual.hpp"
    1214#include "ClpFactorization.hpp"
    1315#include "ClpNonLinearCost.hpp"
    1416#include "ClpLinearObjective.hpp"
    1517#include "ClpConstraint.hpp"
     18#include "ClpPresolve.hpp"
    1619#include "ClpQuadraticObjective.hpp"
    1720#include "CoinPackedMatrix.hpp"
     
    26472650     return returnCode;
    26482651}
     2652// May use a cut approach for solving any LP
     2653int
     2654ClpSimplexNonlinear::primalDualCuts(char * rowsIn, int startUp,int algorithm)
     2655{
     2656  if (!rowsIn) {
     2657    if (algorithm>0)
     2658      return ClpSimplex::primal(startUp);
     2659    else
     2660      //return static_cast<ClpSimplexDual *>(static_cast<ClpSimplex *>(this))->dual(startUp);
     2661      return ClpSimplex::dual(startUp);
     2662  } else {
     2663    int numberUsed=0;
     2664    int rowsThreshold=CoinMax(100,numberRows_/2);
     2665    //int rowsTry=CoinMax(50,numberRows_/4);
     2666    // Just add this number of rows each time in small problem
     2667    int smallNumberRows = 2 * numberColumns_;
     2668    smallNumberRows = CoinMin(smallNumberRows, numberRows_ / 20);
     2669    // We will need arrays to choose rows to add
     2670    double * weight = new double [numberRows_];
     2671    int * sort = new int [numberRows_+numberColumns_];
     2672    int * whichColumns = sort+numberRows_;
     2673    int numberSort = 0;
     2674    for (int i=0;i<numberRows_;i++) {
     2675      if (rowsIn[i])
     2676        numberUsed++;
     2677    }
     2678    if (numberUsed) {
     2679      // normal
     2680    } else {
     2681      // If non slack use that information
     2682      int numberBinding=0;
     2683      numberPrimalInfeasibilities_=0;
     2684      sumPrimalInfeasibilities_=0.0;
     2685      memset(rowActivity_, 0, numberRows_ * sizeof(double));
     2686      times(1.0, columnActivity_, rowActivity_);
     2687      for (int i=0;i<numberRows_;i++) {
     2688        double lowerDifference = rowActivity_[i]-rowLower_[i];
     2689        double upperDifference = rowActivity_[i]-rowUpper_[i];
     2690        if (lowerDifference<-10*primalTolerance_||upperDifference>10.0*primalTolerance_) {
     2691          numberPrimalInfeasibilities_++;
     2692          if (lowerDifference<0.0)
     2693            sumPrimalInfeasibilities_ -= lowerDifference;
     2694          else
     2695            sumPrimalInfeasibilities_ += upperDifference;
     2696          rowsIn[i]=1;
     2697        } else if (getRowStatus(i)!=basic) {
     2698          numberBinding++;
     2699          rowsIn[i]=1;
     2700        }
     2701      }
     2702      if (numberBinding<rowsThreshold) {
     2703        // try
     2704      } else {
     2705        // random?? - or give up
     2706        // Set up initial list
     2707        numberSort = 0;
     2708        for (int iRow = 0; iRow < numberRows_; iRow++) {
     2709          weight[iRow] = 1.123e50;
     2710          if (rowLower_[iRow] == rowUpper_[iRow]) {
     2711            sort[numberSort++] = iRow;
     2712            weight[iRow] = 0.0;
     2713          }
     2714        }
     2715        numberSort /= 2;
     2716        // and pad out with random rows
     2717        double ratio = ((double)(smallNumberRows - numberSort)) / ((double) numberRows_);
     2718        for (int iRow = 0; iRow < numberRows_; iRow++) {
     2719          if (weight[iRow] == 1.123e50 && CoinDrand48() < ratio)
     2720            sort[numberSort++] = iRow;
     2721        }
     2722        // sort
     2723        CoinSort_2(weight, weight + numberRows_, sort);
     2724        numberSort = CoinMin(numberRows_, smallNumberRows);
     2725        memset(rowsIn, 0, numberRows_);
     2726        for (int iRow = 0; iRow < numberSort; iRow++)
     2727          rowsIn[sort[iRow]] = 1;
     2728      }
     2729    }
     2730    // see if feasible
     2731    numberPrimalInfeasibilities_=0;
     2732    memset(rowActivity_, 0, numberRows_ * sizeof(double));
     2733    times(1.0, columnActivity_, rowActivity_);
     2734    for (int i=0;i<numberRows_;i++) {
     2735      double lowerDifference = rowActivity_[i]-rowLower_[i];
     2736      double upperDifference = rowActivity_[i]-rowUpper_[i];
     2737      if (lowerDifference<-10*primalTolerance_||upperDifference>10.0*primalTolerance_) {
     2738        if (lowerDifference<0.0)
     2739          sumPrimalInfeasibilities_ -= lowerDifference;
     2740        else
     2741          sumPrimalInfeasibilities_ += upperDifference;
     2742        numberPrimalInfeasibilities_++;
     2743      }
     2744    }
     2745    printf("Initial infeasibilities - %g (%d)\n",
     2746           sumPrimalInfeasibilities_,numberPrimalInfeasibilities_);
     2747    // Just do this number of passes
     2748    int maxPass = 50;
     2749    // And take out slack rows until this pass
     2750    int takeOutPass = 30;
     2751    int iPass;
     2752   
     2753    const int * start = this->clpMatrix()->getVectorStarts();
     2754    const int * length = this->clpMatrix()->getVectorLengths();
     2755    const int * row = this->clpMatrix()->getIndices();
     2756    problemStatus_=1;
     2757    for (iPass = 0; iPass < maxPass; iPass++) {
     2758      printf("Start of pass %d\n", iPass);
     2759      int numberSmallColumns = 0;
     2760      for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
     2761        int n = 0;
     2762        for (CoinBigIndex j = start[iColumn]; j < start[iColumn] + length[iColumn]; j++) {
     2763          int iRow = row[j];
     2764          if (rowsIn[iRow])
     2765            n++;
     2766        }
     2767        if (n)
     2768          whichColumns[numberSmallColumns++] = iColumn;
     2769      }
     2770      numberSort=0;
     2771      for (int i=0;i<numberRows_;i++) {
     2772        if (rowsIn[i])
     2773          sort[numberSort++]=i;
     2774      }
     2775      // Create small problem
     2776      ClpSimplex small(this, numberSort, sort, numberSmallColumns, whichColumns);
     2777      printf("Small model has %d rows, %d columns and %d elements\n",
     2778             small.numberRows(),small.numberColumns(),small.getNumElements());
     2779      small.setFactorizationFrequency(100 + numberSort / 200);
     2780      // Solve
     2781      small.setLogLevel(CoinMax(0,logLevel()-1));
     2782      if (iPass>20) {
     2783        if (sumPrimalInfeasibilities_>1.0e-1) {
     2784          small.dual();
     2785        } else {
     2786          small.primal(1);
     2787        }
     2788      } else {
     2789        // presolve
     2790#if 0
     2791        ClpSolve::SolveType method;
     2792        ClpSolve::PresolveType presolveType = ClpSolve::presolveOn;
     2793        ClpSolve solveOptions;
     2794        solveOptions.setPresolveType(presolveType, 5);
     2795        if (sumPrimalInfeasibilities_>1.0e-1)
     2796          method = ClpSolve::useDual;
     2797        else
     2798          method = ClpSolve::usePrimalorSprint;
     2799        solveOptions.setSolveType(method);
     2800        small.initialSolve(solveOptions);
     2801#else
     2802#if 1
     2803        ClpPresolve * pinfo = new ClpPresolve();
     2804        ClpSimplex * small2 = pinfo->presolvedModel(small,1.0e-5);
     2805        assert (small2);
     2806        if (sumPrimalInfeasibilities_>1.0e-1) {
     2807          small2->dual();
     2808        } else {
     2809          small2->primal(1);
     2810        }
     2811        pinfo->postsolve(true);
     2812        delete pinfo;
     2813#else
     2814        char * types = new char[small.numberRows()+small.numberColumns()];
     2815        memset(types,0,small.numberRows()+small.numberColumns());
     2816        if (sumPrimalInfeasibilities_>1.0e-1) {
     2817          small.miniSolve(types,types+small.numberRows(),
     2818                          -1,0);
     2819        } else {
     2820          small.miniSolve(types,types+small.numberRows(),
     2821                          1,1);
     2822        }
     2823        delete [] types;
     2824#endif
     2825        if (small.sumPrimalInfeasibilities()>1.0)
     2826          small.primal(1);
     2827#endif
     2828      }
     2829      bool dualInfeasible = (small.status() == 2);
     2830      // move solution back
     2831      const double * smallSolution = small.primalColumnSolution();
     2832      for (int j = 0; j < numberSmallColumns; j++) {
     2833        int iColumn = whichColumns[j];
     2834        columnActivity_[iColumn] = smallSolution[j];
     2835        this->setColumnStatus(iColumn, small.getColumnStatus(j));
     2836      }
     2837      for (int iRow = 0; iRow < numberSort; iRow++) {
     2838        int kRow = sort[iRow];
     2839        this->setRowStatus(kRow, small.getRowStatus(iRow));
     2840      }
     2841      // compute full solution
     2842      memset(rowActivity_, 0, numberRows_ * sizeof(double));
     2843      times(1.0, columnActivity_, rowActivity_);
     2844      if (iPass != maxPass - 1) {
     2845        // Mark row as not looked at
     2846        for (int iRow = 0; iRow < numberRows_; iRow++)
     2847          weight[iRow] = 1.123e50;
     2848        // Look at rows already in small problem
     2849        int iSort;
     2850        int numberDropped = 0;
     2851        int numberKept = 0;
     2852        int numberBinding = 0;
     2853        numberPrimalInfeasibilities_ = 0;
     2854        sumPrimalInfeasibilities_ = 0.0;
     2855        bool allFeasible = small.numberIterations()==0;
     2856        for (iSort = 0; iSort < numberSort; iSort++) {
     2857          int iRow = sort[iSort];
     2858          //printf("%d %g %g\n",iRow,rowActivity_[iRow],small.primalRowSolution()[iSort]);
     2859          if (getRowStatus(iRow) == ClpSimplex::basic) {
     2860            // Basic - we can get rid of if early on
     2861            if (iPass < takeOutPass && !dualInfeasible) {
     2862              double infeasibility = CoinMax(rowActivity_[iRow] - rowUpper_[iRow],
     2863                                             rowLower_[iRow] - rowActivity_[iRow]);
     2864              weight[iRow] = -infeasibility;
     2865              if (infeasibility > primalTolerance_&&!allFeasible) {
     2866                numberPrimalInfeasibilities_++;
     2867                sumPrimalInfeasibilities_ += infeasibility;
     2868              } else {
     2869                weight[iRow] = 1.0;
     2870                numberDropped++;
     2871              }
     2872            } else {
     2873              // keep
     2874              weight[iRow] = -1.0e40;
     2875              numberKept++;
     2876            }
     2877          } else {
     2878            // keep
     2879            weight[iRow] = -1.0e50;
     2880            numberKept++;
     2881            numberBinding++;
     2882          }
     2883        }
     2884        // Now rest
     2885        for (int iRow = 0; iRow < numberRows_; iRow++) {
     2886          sort[iRow] = iRow;
     2887          if (weight[iRow] == 1.123e50) {
     2888            // not looked at yet
     2889            double infeasibility = CoinMax(rowActivity_[iRow] - rowUpper_[iRow],
     2890                                           rowLower_[iRow] - rowActivity_[iRow]);
     2891            weight[iRow] = -infeasibility;
     2892            if (infeasibility > primalTolerance_) {
     2893              numberPrimalInfeasibilities_++;
     2894              sumPrimalInfeasibilities_ += infeasibility;
     2895            }
     2896          }
     2897        }
     2898        // sort
     2899        CoinSort_2(weight, weight + numberRows_, sort);
     2900        numberSort = CoinMin(numberRows_, smallNumberRows + numberKept);
     2901        memset(rowsIn, 0, numberRows_);
     2902        for (int iRow = 0; iRow < numberSort; iRow++)
     2903          rowsIn[sort[iRow]] = 1;
     2904        printf("%d rows binding, %d rows kept, %d rows dropped - new size %d rows, %d columns\n",
     2905               numberBinding, numberKept, numberDropped, numberSort, numberSmallColumns);
     2906        printf("%d rows are infeasible - sum is %g\n", numberPrimalInfeasibilities_,
     2907               sumPrimalInfeasibilities_);
     2908        if (!numberPrimalInfeasibilities_) {
     2909          problemStatus_=0;
     2910          printf("Exiting as looks optimal\n");
     2911          break;
     2912        }
     2913        numberPrimalInfeasibilities_ = 0;
     2914        sumPrimalInfeasibilities_ = 0.0;
     2915        for (iSort = 0; iSort < numberSort; iSort++) {
     2916          if (weight[iSort] > -1.0e30 && weight[iSort] < -1.0e-8) {
     2917            numberPrimalInfeasibilities_++;
     2918            sumPrimalInfeasibilities_ += -weight[iSort];
     2919          }
     2920        }
     2921        printf("in small model %d rows are infeasible - sum is %g\n", numberPrimalInfeasibilities_,
     2922               sumPrimalInfeasibilities_);
     2923      } else {
     2924        // out of passes
     2925        problemStatus_=-1;
     2926      }
     2927    }
     2928    delete [] weight;
     2929    delete [] sort;
     2930    return 0;
     2931  }
     2932}
    26492933// A sequential LP method
    26502934int
    2651 ClpSimplexNonlinear::primalSLP(int numberPasses, double deltaTolerance)
     2935ClpSimplexNonlinear::primalSLP(int numberPasses, double deltaTolerance,
     2936                int otherOptions)
    26522937{
    26532938     // Are we minimizing or maximizing
     
    26902975     objective_ = new ClpLinearObjective(NULL, numberColumns);
    26912976     double * objective = this->objective();
    2692 
     2977     // See if user wants to use cuts
     2978     char * rowsIn=NULL;
     2979     if ((otherOptions&1)!=0||numberPasses<0) {
     2980       numberPasses=abs(numberPasses);
     2981       rowsIn=new char[numberRows_];
     2982       memset(rowsIn,0,numberRows_);
     2983     }
    26932984     // get feasible
    26942985     if (this->status() < 0 || numberPrimalInfeasibilities())
    2695           ClpSimplex::primal(1);
     2986       primalDualCuts(rowsIn,1,1);
    26962987     // still infeasible
    2697      if (numberPrimalInfeasibilities()) {
     2988     if (problemStatus_) {
    26982989          delete [] listNonLinearColumn;
    26992990          return 0;
     
    32653556          if (saveLogLevel == 1)
    32663557               setLogLevel(0);
    3267           ClpSimplex::primal(1);
     3558          primalDualCuts(rowsIn,1,1);
    32683559          algorithm_ = 1;
    32693560          setLogLevel(saveLogLevel);
     
    33393630CoinMemcpyN(trueObjective->gradient(this, solution, offset, true, 2),   numberColumns,
    33403631            objective);
    3341 ClpSimplex::primal(1);
     3632primalDualCuts(rowsIn,1,1);
    33423633delete objective_;
     3634delete [] rowsIn;
    33433635objective_ = trueObjective;
    33443636// redo values
  • trunk/Clp/src/ClpSimplexNonlinear.hpp

    r1665 r2025  
    4141         Using a semi-trust region approach as for pooling problem
    4242         This is in because I have it lying around
    43 
    4443     */
    45      int primalSLP(int numberPasses, double deltaTolerance);
     44     int primalSLP(int numberPasses, double deltaTolerance,
     45                int otherOptions=0);
     46     /// May use a cut approach for solving any LP
     47     int primalDualCuts(char * rowsIn, int startUp, int algorithm);
    4648     /** Primal algorithm for nonlinear constraints
    4749         Using a semi-trust region approach as for pooling problem
  • trunk/Clp/src/ClpSimplexOther.cpp

    r1972 r2025  
    714714                    else
    715715                         fprintf(fp, " UL C%7.7d", iColumn);
     716                    // Dummy row name if values
     717                    if (writeValues)
     718                         fprintf(fp, "      _dummy_");
     719               } else if( (getColumnStatus(iColumn) == ClpSimplex::superBasic||
     720                           getColumnStatus(iColumn) == ClpSimplex::isFree)&&
     721                          writeValues) {
     722                    printit = true;
     723                    if (lengthNames_)
     724                         fprintf(fp, " BS %s", columnNames_[iColumn].c_str());
     725                    else
     726                         fprintf(fp, " BS C%7.7d", iColumn);
    716727                    // Dummy row name if values
    717728                    if (writeValues)
     
    1044410455  delete [] rowLowerX;
    1044510456}
    10446 // Scale real objective
     10457// mini presolve and solve
    1044710458void
    10448 ClpSimplex::scaleRealObjective(double multiplier)
     10459ClpSimplex::miniSolve(char * rowType, char *columnType,int algorithm, int startUp)
    1044910460{
    10450   double * obj = objective();
    10451   for (int i=0;i<numberColumns_;i++)
    10452     obj[i] *= multiplier;
    10453   setObjectiveOffset(multiplier*objectiveOffset());
     10461  listInfo * info = NULL;
     10462  ClpSimplex * small2 = miniPresolve(rowType,columnType,reinterpret_cast<void **>(&info));
     10463  if (algorithm<0)
     10464    small2->dual(startUp);
     10465  else
     10466    small2->primal(startUp);
     10467  miniPostsolve(small2,info);
     10468  delete info;
    1045410469}
    10455 // Scale so no RHS (abs not infinite) > value
    10456 void
    10457 ClpSimplex::scaleRealRhs(double maxValue, double killIfSmaller)
    10458 {
    10459   ClpPackedMatrix * matrix = dynamic_cast<ClpPackedMatrix *>(matrix_);
    10460   if (matrix&&maxValue>1.0) {
    10461     matrix->allElementsInRange(this,
    10462                                1.0e-12,1.0e12,8);
    10463     double * element = matrix->getMutableElements();
    10464     const int * row = matrix->getIndices();
    10465     const int * columnLength = matrix->getVectorLengths();
    10466     const CoinBigIndex * columnStart = matrix->getVectorStarts();
    10467     for (int iRow=0;iRow<numberRows_;iRow++) {
    10468       double lowerValue=rowLower_[iRow];
    10469       if (lowerValue<-1.0e15)
    10470         lowerValue=-1.0;
    10471       double upperValue=rowUpper_[iRow];
    10472       if (upperValue>1.0e15)
    10473         upperValue=1.0;
    10474       double scale = CoinMax(fabs(lowerValue),fabs(upperValue));
    10475       if (scale>maxValue) {
    10476         dual_[iRow]=maxValue/scale;
    10477       } else {
    10478         dual_[iRow]=1.0;
    10479       }
    10480       if (lowerValue>=-1.0e15)
    10481         rowLower_[iRow] *= dual_[iRow];
    10482       if (upperValue<=1.0e15)
    10483         rowUpper_[iRow] *= dual_[iRow];
    10484     }
    10485     for (int iColumn=0;iColumn<numberColumns_;iColumn++) {
    10486       for (CoinBigIndex j = columnStart[iColumn];
    10487            j < columnStart[iColumn] + columnLength[iColumn]; j++) {
    10488         int iRow = row[j];
    10489         element[j] *= dual_[iRow];
    10490       }
    10491     }
    10492     matrix->allElementsInRange(this,
    10493                                killIfSmaller,1.0e22,15);
    10494   } else {
    10495     printf("Can't scale rhs as not ClpPackedMatrix\n");
    10496   }
    10497 }
  • trunk/Clp/src/ClpSolve.cpp

    r2015 r2025  
    4444#include "AbcSimplexFactorization.hpp"
    4545#endif
     46#include "CoinStructuredModel.hpp"
    4647double zz_slack_value=0.0;
    4748
     
    540541#endif
    541542#ifdef ABC_INHERIT
    542 void
     543AbcSimplex *
    543544ClpSimplex::dealWithAbc(int solveType, int startUp,
    544545                        bool interrupt)
    545546{
     547  bool keepAbc = startUp>=10;
     548  startUp=startUp%10;
     549  AbcSimplex * abcModel2 = NULL;
    546550  if (!this->abcState()||!numberRows_||!numberColumns_) {
     551    //this->readBasis("aaa.bas");
    547552    if (!solveType)
    548553      this->dual(0);
    549     else
     554    else if (solveType==1)
    550555      this->primal(startUp ? 1 : 0);
     556    //this->writeBasis("a.bas",true);
    551557  } else {
    552     AbcSimplex * abcModel2=new AbcSimplex(*this);
     558    abcModel2=new AbcSimplex(*this);
    553559    if (interrupt)
    554560      currentAbcModel = abcModel2;
     
    559565    //abcModel2->startPermanentArrays();
    560566    int crashState=abcModel2->abcState()&(256+512+1024);
    561     abcModel2->setAbcState(CLP_ABC_WANTED|crashState);
     567    abcModel2->setAbcState(CLP_ABC_WANTED|crashState|(abcModel2->abcState()&15));
    562568    int ifValuesPass=startUp ? 1 : 0;
    563569    // temp
     
    595601        abcModel2->setupDualValuesPass(this->dualRowSolution(),
    596602                                       this->primalColumnSolution(),type);
    597       } else {
     603      } else if (solveType==1) {
    598604        ifValuesPass=1;
    599605        abcModel2->setStateOfProblem(abcModel2->stateOfProblem() | VALUES_PASS);
     
    674680    if (!solveType) {
    675681      abcModel2->ClpSimplex::doAbcDual();
    676     } else {
     682    } else if (solveType==1) {
    677683      int saveOptions=abcModel2->specialOptions();
    678684      if (startUp==2)
     
    682688    }
    683689#if ABC_INSTRUMENT
    684     double timeCpu=CoinCpuTime()-startTimeCpu;
    685     double timeElapsed=CoinGetTimeOfDay()-startTimeElapsed;
    686     sprintf(line,"Cpu time for %s (%d rows, %d columns %d elements) %g elapsed %g ratio %g - %d iterations",
    687             this->problemName().c_str(),this->numberRows(),this->numberColumns(),
    688             this->getNumElements(),
    689             timeCpu,timeElapsed,timeElapsed ? timeCpu/timeElapsed : 1.0,abcModel2->numberIterations());
    690     handler_->message(CLP_GENERAL, messages_)
    691       << line
    692       << CoinMessageEol;
     690    if (solveType<2) {
     691      double timeCpu=CoinCpuTime()-startTimeCpu;
     692      double timeElapsed=CoinGetTimeOfDay()-startTimeElapsed;
     693      sprintf(line,"Cpu time for %s (%d rows, %d columns %d elements) %g elapsed %g ratio %g - %d iterations",
     694              this->problemName().c_str(),this->numberRows(),this->numberColumns(),
     695              this->getNumElements(),
     696              timeCpu,timeElapsed,timeElapsed ? timeCpu/timeElapsed : 1.0,abcModel2->numberIterations());
     697      handler_->message(CLP_GENERAL, messages_)
     698        << line
     699        << CoinMessageEol;
    693700#if ABC_INSTRUMENT>1
    694     {
    695       int n;
    696       n=0;
    697       for (int i=0;i<20;i++)
    698         n+= abcPricing[i];
    699       printf("CCSparse pricing done %d times",n);
    700       int n2=0;
    701       for (int i=0;i<20;i++)
    702         n2+= abcPricingDense[i];
    703       if (n2)
    704         printf(" and dense pricing done %d times\n",n2);
    705       else
    706         printf("\n");
    707       n=0;
    708       printf ("CCS");
    709       for (int i=0;i<19;i++) {
    710         if (abcPricing[i]) {
     701      {
     702        int n;
     703        n=0;
     704        for (int i=0;i<20;i++)
     705          n+= abcPricing[i];
     706        printf("CCSparse pricing done %d times",n);
     707        int n2=0;
     708        for (int i=0;i<20;i++)
     709          n2+= abcPricingDense[i];
     710        if (n2)
     711          printf(" and dense pricing done %d times\n",n2);
     712        else
     713          printf("\n");
     714        n=0;
     715        printf ("CCS");
     716        for (int i=0;i<19;i++) {
     717          if (abcPricing[i]) {
     718            if (n==5) {
     719              n=0;
     720              printf("\nCCS");
     721            }
     722            n++;
     723            printf("(%d els,%d times) ",i,abcPricing[i]);
     724          }
     725        }
     726        if (abcPricing[19]) {
    711727          if (n==5) {
    712728            n=0;
     
    714730          }
    715731          n++;
    716           printf("(%d els,%d times) ",i,abcPricing[i]);
     732          printf("(>=19 els,%d times) ",abcPricing[19]);
    717733        }
    718       }
    719       if (abcPricing[19]) {
    720         if (n==5) {
    721           n=0;
    722           printf("\nCCS");
    723         }
    724         n++;
    725         printf("(>=19 els,%d times) ",abcPricing[19]);
    726       }
    727       if (n2) {
    728         printf ("CCD");
    729         for (int i=0;i<19;i++) {
    730           if (abcPricingDense[i]) {
    731             if (n==5) {
    732               n=0;
    733               printf("\nCCD");
     734        if (n2) {
     735          printf ("CCD");
     736          for (int i=0;i<19;i++) {
     737            if (abcPricingDense[i]) {
     738              if (n==5) {
     739                n=0;
     740                printf("\nCCD");
     741              }
     742              n++;
     743              int k1=(numberRows_/16)*i;;
     744              int k2=CoinMin(numberRows_,k1+(numberRows_/16)-1);
     745              printf("(%d-%d els,%d times) ",k1,k2,abcPricingDense[i]);
    734746            }
    735             n++;
    736             int k1=(numberRows_/16)*i;;
    737             int k2=CoinMin(numberRows_,k1+(numberRows_/16)-1);
    738             printf("(%d-%d els,%d times) ",k1,k2,abcPricingDense[i]);
    739747          }
    740748        }
     749        printf("\n");
    741750      }
    742       printf("\n");
     751      instrument_print();
     752#endif
    743753    }
    744     instrument_print();
    745 #endif
    746754#endif
    747755    abcModel2->moveStatusToClp(this);
     756#if 1
     757    if (!problemStatus_) {
     758      double offset;
     759      CoinMemcpyN(this->objectiveAsObject()->gradient(this,
     760                                                      this->primalColumnSolution(), offset, true),
     761                  numberColumns_, this->dualColumnSolution());
     762      this->clpMatrix()->transposeTimes(-1.0,
     763                                        this->dualRowSolution(),
     764                                        this->dualColumnSolution());
     765      memset(this->primalRowSolution(), 0, numberRows_ * sizeof(double));
     766      this->clpMatrix()->times(1.0,
     767                               this->primalColumnSolution(),
     768                               this->primalRowSolution());
     769      this->checkSolutionInternal();
     770      if (sumDualInfeasibilities_>100.0*dualTolerance_) {
     771        printf("internal check on duals failed %d %g\n",
     772               numberDualInfeasibilities_,sumDualInfeasibilities_);
     773      } else {
     774        sumDualInfeasibilities_=0.0;
     775        numberDualInfeasibilities_=0;
     776      }
     777      if (sumPrimalInfeasibilities_>100.0*primalTolerance_) {
     778        printf("internal check on primals failed %d %g\n",
     779               numberPrimalInfeasibilities_,sumPrimalInfeasibilities_);
     780      } else {
     781        sumPrimalInfeasibilities_=0.0;
     782        numberPrimalInfeasibilities_=0;
     783      }
     784      problemStatus_=0;
     785    }
     786#endif
    748787    //ClpModel::stopPermanentArrays();
    749788    this->setSpecialOptions(this->specialOptions()&~65536);
     
    760799    moreSpecialOptions_ &= ~16384;
    761800    //this->setNumberIterations(abcModel2->numberIterations()+this->numberIterations());
    762     delete abcModel2;
     801    if (!keepAbc) {
     802      delete abcModel2;
     803      abcModel2=NULL;
     804    }
    763805  }
     806  return abcModel2;
    764807}
    765808#endif
     
    785828     double time1 = CoinCpuTime();
    786829     double timeX = time1;
    787      double time2;
     830     double time2=0.0;
    788831     ClpMatrixBase * saveMatrix = NULL;
    789832     ClpObjective * savedObjective = NULL;
     
    9811024     }
    9821025     if (method != ClpSolve::useDual && method != ClpSolve::useBarrier
     1026         && method != ClpSolve::tryBenders && method != ClpSolve::tryDantzigWolfe
    9831027               && method != ClpSolve::useBarrierNoCross) {
    9841028          switch (options.getSpecialOption(1)) {
     
    10411085               doCrash = 0;
    10421086               doSprint = 0;
    1043                if (options.getExtraInfo(1) > 0)
     1087               if (options.getExtraInfo(1) )
    10441088                    doSlp = options.getExtraInfo(1);
    10451089               break;
     
    10531097               abort();
    10541098          }
    1055      } else {
     1099     } else if (method != ClpSolve::tryBenders && method != ClpSolve::tryDantzigWolfe) {
    10561100          // Dual
    10571101          switch (options.getSpecialOption(0)) {
     
    10781122               abort();
    10791123          }
     1124     } else {
     1125       // decomposition
    10801126     }
    10811127#ifndef NO_RTTI
     
    10891135     if (quadraticObj) {
    10901136          doSprint = 0;
    1091           doIdiot = 0;
     1137          //doIdiot = 0;
    10921138          // off
    10931139          if (method == ClpSolve::useBarrier)
     
    11291175          if (model2->getRowStatus(iRow) == basic)
    11301176               numberRowsBasic++;
    1131      if (numberRowsBasic < numberRows) {
     1177     if (numberRowsBasic < numberRows && objective_->type()<2) {
    11321178          doIdiot = 0;
    11331179          doCrash = 0;
     
    12221268                    const double * objective = model2->objective();
    12231269                    int nObj=0;
     1270                    int nFree=0;
    12241271                    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
    12251272                         if (columnLower[iColumn] < -1.0e10 && columnUpper[iColumn] > 1.0e10) {
    1226                               doSprint = 0;
    1227                               doIdiot = 0;
    1228                               break;
     1273                           nFree++;
    12291274                         } else if (objective[iColumn]) {
    12301275                           nObj++;
     
    12321277                    }
    12331278                    if (nObj*10<numberColumns)
     1279                      doSprint=0;
     1280                    if (nFree)
     1281                      doIdiot=0;
     1282                    if (nFree*10>numberRows)
    12341283                      doSprint=0;
    12351284                    int nPasses = 0;
     
    13251374          }
    13261375     }
     1376     if (method == ClpSolve::tryBenders) {
     1377       // Now build model
     1378       int lengthNames=model2->lengthNames();
     1379       model2->setLengthNames(0);
     1380       CoinModel * build = model2->createCoinModel();
     1381       model2->setLengthNames(lengthNames);
     1382       CoinStructuredModel benders;
     1383       build->convertMatrix();
     1384       int numberBlocks = options.independentOption(0);
     1385       benders.setMessageHandler(handler_);
     1386       numberBlocks=benders.decompose(*build,2,numberBlocks,NULL);
     1387       delete build;
     1388       //exit(0);
     1389       if (numberBlocks) {
     1390         options.setIndependentOption(1,1); // don't do final clean up
     1391         model2->solveBenders(&benders,options);
     1392         //move solution
     1393         method=ClpSolve::notImplemented;
     1394         time2 = CoinCpuTime();
     1395         timeCore = time2 - timeX;
     1396         handler_->message(CLP_INTERVAL_TIMING, messages_)
     1397           << "Crossover" << timeCore << time2 - time1
     1398           << CoinMessageEol;
     1399         timeX = time2;
     1400       } else {
     1401         printf("No structure\n");
     1402         method=ClpSolve::useDual;
     1403       }
     1404     } else if (method == ClpSolve::tryDantzigWolfe) {
     1405         abort();
     1406     }
    13271407     if (method == ClpSolve::usePrimalorSprint) {
    13281408          if (doSprint < 0) {
     
    17991879                    // Allow for scaling
    18001880                    info.setStrategy(32 | info.getStrategy());
    1801                     info.crash(nPasses, model2->messageHandler(), model2->messagesPointer());
     1881                    int saveScalingFlag=model2->scalingFlag();
     1882                    if (objective_->type()==2)
     1883                      model2->scaling(0);
     1884                    info.crash(nPasses, model2->messageHandler(), model2->messagesPointer(),(objective_->type() <2));
     1885                      model2->scaling(saveScalingFlag);
    18021886#endif
    18031887                    time2 = CoinCpuTime();
     
    18921976          }
    18931977#ifndef SLIM_CLP
    1894           if (doSlp > 0 && objective_->type() == 2) {
     1978          if (doSlp && objective_->type() == 2) {
    18951979               model2->nonlinearSLP(doSlp, 1.0e-5);
    18961980          }
     
    21772261                    sort[i] = i + originalNumberColumns;
    21782262          }
     2263          // put in free
    21792264          // maybe a solution there already
    21802265          for (iColumn = 0; iColumn < originalNumberColumns; iColumn++) {
    2181                if (model2->getColumnStatus(iColumn) == basic)
     2266               if (model2->getColumnStatus(iColumn) == basic||
     2267                   (columnLower[iColumn]<-1.0e30&&columnUpper[iColumn]>1.0e30))
    21822268                    sort[numberSort++] = iColumn;
    21832269          }
     
    22812367                         small.dealWithAbc(0,0);
    22822368#else
    2283                       if (iPass||!numberArtificials) 
     2369                      if (iPass||!numberArtificials)
    22842370                         small.primal(1);
    22852371                      else
     
    25472633#ifdef COIN_HAS_WSMP
    25482634          case 2: {
    2549                ClpCholeskyWssmp * cholesky = new ClpCholeskyWssmp(CoinMax(100, model2->numberRows() / 10));
    2550                barrier.setCholesky(cholesky);
    2551                assert (!doKKT);
     2635               if (!doKKT) {
     2636                    ClpCholeskyWssmp * cholesky = new ClpCholeskyWssmp(CoinMax(100, model2->numberRows() / 10));
     2637                    barrier.setCholesky(cholesky);
     2638               } else {
     2639                    ClpCholeskyWssmpKKT * cholesky = new ClpCholeskyWssmpKKT(CoinMax(100, model2->numberRows() / 10));
     2640                    barrier.setCholesky(cholesky);
     2641               }
    25522642          }
    25532643          break;
     
    25842674#ifdef COIN_HAS_MUMPS
    25852675          case 6: {
    2586                ClpCholeskyMumps * cholesky = new ClpCholeskyMumps();
    2587                barrier.setCholesky(cholesky);
    2588                assert (!doKKT);
     2676               if (!doKKT) {
     2677                    ClpCholeskyMumps * cholesky = new ClpCholeskyMumps();
     2678                    barrier.setCholesky(cholesky);
     2679               } else {
     2680                    ClpCholeskyMumps * cholesky = new ClpCholeskyMumps();
     2681                    cholesky->setKKT(true);
     2682                    barrier.setCholesky(cholesky);
     2683               }
    25892684          }
    25902685          break;
     
    29042999                                   saveUpper = NULL;
    29053000                              }
    2906                               int numberRows = model2->numberRows();
    2907                               int numberColumns = model2->numberColumns();
     3001                              //int numberRows = model2->numberRows();
     3002                              //int numberColumns = model2->numberColumns();
    29083003#ifdef ABC_INHERIT
    29093004                              model2->checkSolution(0);
     
    30623157          abort();
    30633158#endif
     3159     } else if (method == ClpSolve::notImplemented) {
     3160         printf("done decomposition\n");
    30643161     } else {
    30653162          assert (method != ClpSolve::automatic); // later
     
    36763773                         // dual - change tolerance
    36773774                         model_->setCurrentDualTolerance(model_->currentDualTolerance() * 1.05);
    3678                          model_->setCurrentPrimalTolerance(model_->currentPrimalTolerance() * 1.05);
    36793775                         // if infeasible increase dual bound
    36803776                         if (model_->dualBound() < 1.0e17) {
     
    38043900     realInfeasibility_[CLP_PROGRESS-1] = value;
    38053901}
     3902// Returns number of primal infeasibilities (if -1) - current if (0)
     3903int
     3904ClpSimplexProgress::numberInfeasibilities(int back) const
     3905{
     3906    return numberInfeasibilities_[CLP_PROGRESS-1-back];
     3907}
    38063908// Modify objective e.g. if dual infeasible in dual
    38073909void
     
    40264128     delete [] blockInfo;
    40274129     // decide what to do
     4130     ClpSolve options;
     4131     options.setIndependentOption(2,100);
    40284132     switch (decomposeType) {
    40294133          // No good
     
    40334137          // DW
    40344138     case 1:
    4035           return solveDW(model);
     4139       return solveDW(model,options);
    40364140          // Benders
    40374141     case 2:
    4038           return solveBenders(model);
     4142       return solveBenders(model,options);
    40394143     }
    40404144     return 0; // to stop compiler warning
     
    40824186     CoinFillN(rowBase, numberRowBlocks, -1);
    40834187     // And row to put it
    4084      int * whichRow = new int [numberRows];
     4188     int * whichRow = new int [numberRows+numberRowBlocks];
    40854189     int * columnBase = new int[numberColumnBlocks];
    40864190     CoinFillN(columnBase, numberColumnBlocks, -1);
    40874191     // And column to put it
    4088      int * whichColumn = new int [numberColumns];
     4192     int * whichColumn = new int [numberColumns+numberColumnBlocks];
    40894193     for (int iBlock = 0; iBlock < numberElementBlocks; iBlock++) {
    40904194          CoinModel * block = coinModel.coinBlock(iBlock);
     
    41024206               rowBase[iRowBlock] = block->numberRows();
    41034207               // Save block number
    4104                whichRow[numberRows-numberRowBlocks+iRowBlock] = iBlock;
     4208               whichRow[numberRows+iRowBlock] = iBlock;
    41054209          } else {
    41064210               assert(rowBase[iRowBlock] == block->numberRows());
     
    41094213               columnBase[iColumnBlock] = block->numberColumns();
    41104214               // Save block number
    4111                whichColumn[numberColumns-numberColumnBlocks+iColumnBlock] = iBlock;
     4215               whichColumn[numberColumns+iColumnBlock] = iBlock;
    41124216          } else {
    41134217               assert(columnBase[iColumnBlock] == block->numberColumns());
     
    41274231          assert (k >= 0);
    41284232          // block number
    4129           int jBlock = whichRow[numberRows-numberRowBlocks+iBlock];
     4233          int jBlock = whichRow[numberRows+iBlock];
    41304234          if (originalOrder) {
    41314235               memcpy(whichRow + n, coinModel.coinBlock(jBlock)->originalRows(), k * sizeof(int));
     
    41434247          if (k) {
    41444248               // block number
    4145                int jBlock = whichColumn[numberColumns-numberColumnBlocks+iBlock];
     4249               int jBlock = whichColumn[numberColumns+iBlock];
    41464250               if (originalOrder) {
    41474251                    memcpy(whichColumn + n, coinModel.coinBlock(jBlock)->originalColumns(),
     
    43444448// Solve using Dantzig-Wolfe decomposition and maybe in parallel
    43454449int
    4346 ClpSimplex::solveDW(CoinStructuredModel * model)
     4450ClpSimplex::solveDW(CoinStructuredModel * model,ClpSolve & options)
    43474451{
    43484452     double time1 = CoinCpuTime();
     
    44414545     // Set offset
    44424546     master.setObjectiveOffset(model->objectiveOffset());
     4547     bool reducePrint = logLevel()==7;
     4548     if (reducePrint) {
     4549       // special
     4550       setLogLevel(1);
     4551       master.setLogLevel(1);
     4552     }
    44434553     kBlock = 0;
    44444554     int masterBlock = -1;
     
    45604670          lastArtificial = master.numberColumns();
    45614671     }
    4562      int maxPass = 500;
     4672     int maxPass = options.independentOption(2);
     4673     if (maxPass<2)
     4674       maxPass=100;
    45634675     int iPass;
    45644676     double lastObjective = 1.0e31;
     
    45994711                            columnAdd, rowAdd, elementAdd);
    46004712     }
     4713     char generalPrint[200];
    46014714     // and resize matrix to double check clp will be happy
    46024715     //master.matrix()->setDimensions(numberMasterRows+numberBlocks,
    46034716     //                  numberMasterColumns+numberBlocks);
    4604      std::cout << "Time to decompose " << CoinCpuTime() - time1 << " seconds" << std::endl;
     4717     sprintf(generalPrint,"Time to decompose %.2f seconds",CoinCpuTime() - time1);
     4718     handler_->message(CLP_GENERAL, messages_)
     4719       << generalPrint
     4720       << CoinMessageEol;
    46054721     for (iPass = 0; iPass < maxPass; iPass++) {
    4606           printf("Start of pass %d\n", iPass);
     4722          sprintf(generalPrint,"Start of pass %d", iPass);
     4723          handler_->message(CLP_GENERAL, messages_)
     4724            << generalPrint
     4725            << CoinMessageEol;
    46074726          // Solve master - may be infeasible
    46084727          //master.scaling(0);
     
    46354754                    }
    46364755               }
    4637                printf("Sum of artificials before solve is %g\n", sumArtificials);
     4756               sprintf(generalPrint,"Sum of artificials before solve is %g", sumArtificials);
     4757               handler_->message(CLP_GENERAL, messages_)
     4758                 << generalPrint
     4759                 << CoinMessageEol;
    46384760          }
    46394761          // scale objective to be reasonable
     
    46554777                         sum -= value - lower[iRow];
    46564778               }
    4657                printf("suminf %g\n", sum);
     4779               printf("** suminf %g\n", sum);
    46584780               lower = master.columnLower();
    46594781               upper = master.columnUpper();
     
    46664788                         sum -= value - lower[iColumn];
    46674789               }
    4668                printf("suminf %g\n", sum);
     4790               printf("** suminf %g\n", sum);
    46694791          }
    46704792          master.primal(1);
     
    46764798                    sumArtificials += solution[i];
    46774799               }
    4678                printf("Sum of artificials after solve is %g\n", sumArtificials);
     4800               printf("** Sum of artificials after solve is %g\n", sumArtificials);
    46794801          }
    46804802          master.scaleObjective(scaleFactor);
    46814803          int problemStatus = master.status(); // do here as can change (delcols)
     4804          if (problemStatus == 2 && master.numberColumns()) {
     4805            master.primal(1);
     4806            if (problemStatus==2) {
     4807              int numberColumns = master.numberColumns();
     4808              double * lower = master.columnLower();
     4809              double * upper = master.columnUpper();
     4810              for (int i=0;i<numberColumns;i++) {
     4811                lower[i]=CoinMax(lower[i],-1.0e10);
     4812                upper[i]=CoinMin(upper[i],1.0e10);
     4813              }
     4814              master.primal(1);
     4815              assert (problemStatus!=2);
     4816            }
     4817          }
    46824818          if (master.numberIterations() == 0 && iPass)
    46834819               break; // finished
     
    47194855               dual = master.infeasibilityRay();
    47204856               deleteDual = true;
    4721                printf("The sum of infeasibilities is %g\n",
     4857               printf("** The sum of infeasibilities is %g\n",
    47224858                      master.sumPrimalInfeasibilities());
    47234859          } else if (!master.numberColumns()) {
     
    47274863                      0, (numberMasterRows + numberBlocks)*sizeof(double));
    47284864          } else {
     4865               master.writeMps("unbounded.mps");
    47294866               abort();
    47304867          }
     
    47584895               double scaleFactor =
    47594896                    sub[iBlock].scaleObjective((sumArtificials > 1.0e-5) ? -1.0e-4 : -1.0e9);
     4897                   
     4898               if (reducePrint)
     4899                   sub[iBlock].setLogLevel(0);
    47604900               if (iPass) {
    47614901                    sub[iBlock].primal();
     
    48134953                         // if elements large then scale?
    48144954                         //if (largest>1.0e8||smallest<1.0e-8)
    4815                          printf("For subproblem %d smallest - %g, largest %g - dj %g\n",
     4955                         sprintf(generalPrint,"For subproblem %d smallest - %g, largest %g - dj %g",
    48164956                                iBlock, smallest, largest, dj);
     4957                         handler_->message(CLP_GENERAL2, messages_)
     4958                           << generalPrint
     4959                           << CoinMessageEol;
    48174960                         if (dj < -1.0e-6 || !iPass) {
    48184961                              // take
     
    48454988                         // if elements large or small then scale?
    48464989                         //if (largest>1.0e8||smallest<1.0e-8)
    4847                          printf("For subproblem ray %d smallest - %g, largest %g - dj %g\n",
     4990                         sprintf(generalPrint,"For subproblem ray %d smallest - %g, largest %g - dj %g",
    48484991                                iBlock, smallest, largest, dj);
     4992                         handler_->message(CLP_GENERAL2, messages_)
     4993                           << generalPrint
     4994                           << CoinMessageEol;
    48494995                         if (dj < -1.0e-6) {
    48504996                              // take
     
    48665012                                 columnAdd, rowAdd, elementAdd);
    48675013     }
    4868      std::cout << "Time at end of D-W " << CoinCpuTime() - time1 << " seconds" << std::endl;
     5014     sprintf(generalPrint,"Time at end of D-W %.2f seconds",CoinCpuTime() - time1);
     5015     handler_->message(CLP_GENERAL, messages_)
     5016       << generalPrint
     5017       << CoinMessageEol;
    48695018     //master.scaling(0);
    48705019     //master.primal(1);
     
    49765125     //for (int iRow=0;iRow<numberRows;iRow++)
    49775126     //setRowStatus(iRow,ClpSimplex::superBasic);
    4978      std::cout << "Time before cleanup of full model " << CoinCpuTime() - time1 << " seconds" << std::endl;
     5127     sprintf(generalPrint,"Time before cleanup of full model %.2f seconds",CoinCpuTime() - time1);
     5128     handler_->message(CLP_GENERAL, messages_)
     5129       << generalPrint
     5130       << CoinMessageEol;
    49795131     primal(1);
    4980      std::cout << "Total time " << CoinCpuTime() - time1 << " seconds" << std::endl;
     5132     sprintf(generalPrint,"Total time %.2f seconds",CoinCpuTime() - time1);
     5133     handler_->message(CLP_GENERAL, messages_)
     5134       << generalPrint
     5135       << CoinMessageEol;
    49815136     delete [] columnCounts;
    49825137     delete [] sol;
     
    49935148     return 0;
    49945149}
     5150static ClpSimplex * deBound(ClpSimplex * oldModel)
     5151{
     5152  ClpSimplex * model = new ClpSimplex(*oldModel);
     5153  int numberRows = model->numberRows();
     5154  CoinPackedMatrix * matrix = model->matrix();
     5155  const int * row = matrix->getIndices();
     5156  const int * columnLength = matrix->getVectorLengths();
     5157  const CoinBigIndex * columnStart = matrix->getVectorStarts();
     5158  double * elementByColumn = matrix->getMutableElements();
     5159  int numberColumns = model->numberColumns();
     5160  double * rowLower = model->rowLower();
     5161  double * rowUpper = model->rowUpper();
     5162  double * columnLower = model->columnLower();
     5163  double * columnUpper = model->columnUpper();
     5164  double * objective = model->objective();
     5165  double * change = new double[CoinMax(numberRows,numberColumns)+numberColumns];
     5166  int * rowStart = new int [2*numberColumns+1];
     5167  memset(change,0,numberRows*sizeof(double));
     5168  // first swap ones with infinite lower bounds
     5169  for (int iColumn=0;iColumn<numberColumns;iColumn++) {
     5170    if (columnLower[iColumn]==-COIN_DBL_MAX&&
     5171        columnUpper[iColumn]!=COIN_DBL_MAX) {
     5172      for (CoinBigIndex j=columnStart[iColumn];j<columnStart[iColumn]
     5173             +columnLength[iColumn];j++)
     5174        elementByColumn[j] *= -1.0;
     5175      objective[iColumn] *= -1.0;
     5176      columnLower[iColumn] = -columnUpper[iColumn];
     5177      columnUpper[iColumn] = COIN_DBL_MAX;
     5178    }
     5179  }
     5180  // Out nonzero LB's
     5181  for (int iColumn=0;iColumn<numberColumns;iColumn++) {
     5182    if (columnLower[iColumn]) {
     5183      double value=columnLower[iColumn];
     5184      for (CoinBigIndex j=columnStart[iColumn];j<columnStart[iColumn]
     5185             +columnLength[iColumn];j++) {
     5186        int iRow=row[j];
     5187        change[iRow] -= value*elementByColumn[j];
     5188      }
     5189    }
     5190  }
     5191  for (int iRow = 0; iRow < numberRows; iRow++) {
     5192    double value=change[iRow];
     5193    if (rowLower[iRow]>-COIN_DBL_MAX)
     5194      rowLower[iRow]-=value;
     5195    if (rowUpper[iRow]<COIN_DBL_MAX)
     5196      rowUpper[iRow]-=value;
     5197  }
     5198  int nExtra=0;
     5199  int * columnNew = rowStart+numberColumns+1;
     5200  for (int iColumn=0;iColumn<numberColumns;iColumn++) {
     5201    if (columnUpper[iColumn]<COIN_DBL_MAX&&columnUpper[iColumn]) {
     5202      columnNew[nExtra]=iColumn;
     5203      change[nExtra++]=columnUpper[iColumn];
     5204      columnUpper[iColumn]=COIN_DBL_MAX;
     5205    }
     5206  }
     5207  double * elementNew = change+numberColumns;
     5208  for (int i=0;i<nExtra;i++) {
     5209    rowStart[i]=i;
     5210    elementNew[i]=1.0;
     5211  }
     5212  rowStart[nExtra]=nExtra;
     5213  model->addRows(nExtra, NULL, change,
     5214                 rowStart, columnNew, elementNew);
     5215  delete [] rowStart;
     5216  delete [] change;
     5217  return model;
     5218}
     5219#ifdef ABC_INHERIT
     5220static void * abc_parallelManager(void * stuff)
     5221{
     5222  AbcPthreadStuff * driver = reinterpret_cast<AbcPthreadStuff *>(stuff);
     5223  int whichThread=driver->whichThread();
     5224  CoinAbcThreadInfo * threadInfo = driver->threadInfoPointer(whichThread);
     5225  threadInfo->status=-1;
     5226  int * which = threadInfo->stuff;
     5227  pthread_barrier_wait(driver->barrierPointer());
     5228#if 0
     5229  int status=-1;
     5230  while (status!=100)
     5231    status=timedWait(driver,1000,2);
     5232  pthread_cond_signal(driver->conditionPointer(1));
     5233  pthread_mutex_unlock(driver->mutexPointer(1,whichThread));
     5234#endif
     5235  // so now mutex_ is locked
     5236  int whichLocked=0;
     5237  while (true) {
     5238    pthread_mutex_t * mutexPointer = driver->mutexPointer(whichLocked,whichThread);
     5239    // wait
     5240    //printf("Child waiting for %d - status %d %d %d\n",
     5241    //     whichLocked,lockedX[0],lockedX[1],lockedX[2]);
     5242#ifdef DETAIL_THREAD
     5243    printf("thread %d about to lock mutex %d\n",whichThread,whichLocked);
     5244#endif
     5245    pthread_mutex_lock (mutexPointer);
     5246    whichLocked++;
     5247    if (whichLocked==3)
     5248      whichLocked=0;
     5249    int unLock=whichLocked+1;
     5250    if (unLock==3)
     5251      unLock=0;
     5252    //printf("child pointer %x status %d\n",threadInfo,threadInfo->status);
     5253    assert(threadInfo->status>=0);
     5254    if (threadInfo->status==1000)
     5255      pthread_exit(NULL);
     5256    int type=threadInfo->status;
     5257    int & returnCode=which[0];
     5258    int iPass=which[1];
     5259    ClpSimplex * clpSimplex = reinterpret_cast<ClpSimplex *>(threadInfo->extraInfo);
     5260    //CoinIndexedVector * array;
     5261    //double dummy;
     5262    switch(type) {
     5263      // dummy
     5264    case 0:
     5265      break;
     5266    case 1:
     5267      if (!clpSimplex->problemStatus()||!iPass)
     5268        returnCode=clpSimplex->dual();
     5269      else
     5270        returnCode=clpSimplex->primal();
     5271      break;
     5272    case 100:
     5273      // initialization
     5274      break;
     5275    }
     5276    threadInfo->status=-1;
     5277#ifdef DETAIL_THREAD
     5278    printf("thread %d about to unlock mutex %d\n",whichThread,unLock);
     5279#endif
     5280    pthread_mutex_unlock (driver->mutexPointer(unLock,whichThread));
     5281  }
     5282}
     5283AbcPthreadStuff::AbcPthreadStuff(int numberThreads)
     5284{
     5285  numberThreads_=numberThreads;
     5286  if (numberThreads>8)
     5287    numberThreads=1;
     5288  // For waking up thread
     5289  memset(mutex_,0,sizeof(mutex_));
     5290  for (int iThread=0;iThread<numberThreads;iThread++) {
     5291    for (int i=0;i<3;i++) {
     5292      pthread_mutex_init(&mutex_[i+3*iThread], NULL);
     5293      if (i<2)
     5294        pthread_mutex_lock(&mutex_[i+3*iThread]);
     5295    }
     5296    threadInfo_[iThread].status = 100;
     5297  }
     5298  //pthread_barrierattr_t attr;
     5299  pthread_barrier_init(&barrier_, /*&attr*/ NULL, numberThreads+1);
     5300  for (int iThread=0;iThread<numberThreads;iThread++) {
     5301    pthread_create(&abcThread_[iThread], NULL, abc_parallelManager, reinterpret_cast<void *>(this));
     5302  }
     5303  pthread_barrier_wait(&barrier_);
     5304  pthread_barrier_destroy(&barrier_);
     5305  for (int iThread=0;iThread<numberThreads;iThread++) {
     5306    threadInfo_[iThread].status = -1;
     5307    threadInfo_[iThread].stuff[3]=1; // idle
     5308    locked_[iThread]=0;
     5309  }
     5310}
     5311AbcPthreadStuff::~AbcPthreadStuff()
     5312{
     5313  for (int iThread=0;iThread<numberThreads_;iThread++) {
     5314    startParallelTask(1000,iThread);
     5315  }
     5316  for (int iThread=0;iThread<numberThreads_;iThread++) {
     5317    pthread_join(abcThread_[iThread],NULL);
     5318    for (int i=0;i<3;i++) {
     5319      pthread_mutex_destroy (&mutex_[i+3*iThread]);
     5320    }
     5321  }
     5322}
     5323// so thread can find out which one it is
     5324int
     5325AbcPthreadStuff::whichThread() const
     5326{
     5327  pthread_t thisThread=pthread_self();
     5328  int whichThread;
     5329  for (whichThread=0;whichThread<numberThreads_;whichThread++) {
     5330    if (pthread_equal(thisThread,abcThread_[whichThread]))
     5331      break;
     5332  }
     5333  assert (whichThread<NUMBER_THREADS+1);
     5334  return whichThread;
     5335}
     5336void
     5337AbcPthreadStuff::startParallelTask(int type, int iThread, void * info)
     5338{
     5339  /*
     5340    first time 0,1 owned by main 2 by child
     5341    at end of cycle should be 1,2 by main 0 by child then 2,0 by main 1 by child
     5342  */
     5343  threadInfo_[iThread].status=type;
     5344  threadInfo_[iThread].extraInfo=info;
     5345  threadInfo_[iThread].stuff[3]=0; // say not idle
     5346#ifdef DETAIL_THREAD
     5347  printf("main doing thread %d about to unlock mutex %d\n",iThread,locked_[iThread]);
     5348#endif
     5349  pthread_mutex_unlock (&mutex_[locked_[iThread]+3*iThread]);
     5350}
     5351int
     5352AbcPthreadStuff::waitParallelTask(int type,int & iThread, bool allowIdle)
     5353{
     5354  bool finished = false;
     5355  if (allowIdle) {
     5356    for (iThread = 0; iThread < numberThreads_; iThread++) {
     5357      if (threadInfo_[iThread].status < 0&&threadInfo_[iThread].stuff[3]) {
     5358        finished = true;
     5359        break;
     5360      }
     5361    }
     5362    if (finished)
     5363      return 0;
     5364  }
     5365  while (!finished) {
     5366    for (iThread = 0; iThread < numberThreads_; iThread++) {
     5367      if (threadInfo_[iThread].status < 0&&!threadInfo_[iThread].stuff[3]) {
     5368        finished = true;
     5369        break;
     5370      }
     5371    }
     5372    if (!finished) {
     5373      struct timespec time1, time2;
     5374      time1.tv_sec = 0;
     5375      time1.tv_nsec = 100000;
     5376      nanosleep(&time1 , &time2);
     5377    }
     5378  }
     5379  int locked=locked_[iThread]+2;
     5380  if (locked>=3)
     5381    locked-=3;
     5382#ifdef DETAIL_THREAD
     5383  printf("Main do thread %d about to lock mutex %d\n",iThread,locked);
     5384#endif
     5385  pthread_mutex_lock(&mutex_[locked+iThread*3]);
     5386  locked_[iThread]++;
     5387  if (locked_[iThread]==3)
     5388    locked_[iThread]=0;
     5389  threadInfo_[iThread].stuff[3]=1; // say idle
     5390  return threadInfo_[iThread].stuff[0];
     5391}
     5392void
     5393AbcPthreadStuff::waitAllTasks()
     5394{
     5395  int nWait=0;
     5396  for (int iThread=0;iThread<numberThreads_;iThread++) {
     5397    int idle = threadInfo_[iThread].stuff[3];
     5398    if (!idle)
     5399      nWait++;
     5400  }
     5401#ifdef DETAIL_THREAD
     5402  printf("Waiting for %d tasks to finish\n",nWait);
     5403#endif
     5404  for (int iThread=0;iThread<nWait;iThread++) {
     5405    int jThread;
     5406    waitParallelTask(0,jThread,false);
     5407#ifdef DETAIL_THREAD
     5408    printf("finished with thread %d\n",jThread);
     5409#endif
     5410  }
     5411}
     5412#endif
    49955413// Solve using Benders decomposition and maybe in parallel
    49965414int
    4997 ClpSimplex::solveBenders(CoinStructuredModel * model)
     5415ClpSimplex::solveBenders(CoinStructuredModel * model,ClpSolve & options)
    49985416{
    49995417     double time1 = CoinCpuTime();
     5418     //ClpSimplex * xxxx = deBound(this);
     5419     //xxxx->writeMps("nobounds.mps");
     5420     //delete xxxx;
    50005421     //int numberColumns = model->numberColumns();
    50015422     int numberRowBlocks = model->numberRowBlocks();
    50025423     int numberColumnBlocks = model->numberColumnBlocks();
    50035424     int numberElementBlocks = model->numberElementBlocks();
     5425     char generalPrint[200];
    50045426     // We already have top level structure
    50055427     CoinModelBlockInfo * blockInfo = new CoinModelBlockInfo [numberElementBlocks];
     
    50895511     const CoinPackedMatrix ** first = new const CoinPackedMatrix * [numberRowBlocks];
    50905512     ClpSimplex * sub = new ClpSimplex [numberBlocks];
    5091      ClpSimplex master;
     5513     ClpSimplex masterModel;
    50925514     // Set offset
    5093      master.setObjectiveOffset(model->objectiveOffset());
     5515     masterModel.setObjectiveOffset(model->objectiveOffset());
    50945516     kBlock = 0;
     5517#define ADD_ARTIFICIALS
     5518#ifdef ADD_ARTIFICIALS
     5519     int * originalSubColumns = new int [numberBlocks];
     5520     for (iBlock = 0; iBlock < numberBlocks; iBlock++)
     5521       originalSubColumns[iBlock]=9999999;
     5522#endif
    50955523     int masterBlock = -1;
    50965524     for (iBlock = 0; iBlock < numberRowBlocks; iBlock++) {
     
    51365564                         }
    51375565                         sub[kBlock].setPerturbation(50);
     5566#ifdef ADD_ARTIFICIALS
     5567                         originalSubColumns[kBlock]=sub[kBlock].numberColumns();
     5568                         if (intParam_[0]<1000000) {
     5569                           printf("** Adding artificials\n");
     5570                           int nRow = sub[kBlock].numberRows();
     5571                           int * addStarts = new int [2*nRow+1];
     5572                           int * addRow = new int[2*nRow];
     5573                           double * addElement = new double[2*nRow];
     5574                           addStarts[0] = 0;
     5575                           int numberArtificials = 0;
     5576                           double penalty =1.0e5;
     5577                           double * addCost = new double [2*nRow];
     5578                           const double * lower = sub[kBlock].rowLower();
     5579                           const double * upper = sub[kBlock].rowUpper();
     5580                           for (int iRow = 0; iRow < nRow; iRow++) {
     5581                             if (lower[iRow] > -1.0e20) {
     5582                               addRow[numberArtificials] = iRow;
     5583                               addElement[numberArtificials] = 1.0;
     5584                               addCost[numberArtificials] = penalty;
     5585                               numberArtificials++;
     5586                               addStarts[numberArtificials] = numberArtificials;
     5587                             }
     5588                             if (upper[iRow] < 1.0e20) {
     5589                               addRow[numberArtificials] = iRow;
     5590                               addElement[numberArtificials] = -1.0;
     5591                               addCost[numberArtificials] = penalty;
     5592                               numberArtificials++;
     5593                               addStarts[numberArtificials] = numberArtificials;
     5594                             }
     5595                           }
     5596                           if (numberArtificials) {
     5597                             sub[kBlock].addColumns(numberArtificials, NULL, NULL, addCost,
     5598                                                    addStarts, addRow, addElement);
     5599                           }
     5600                           delete [] addStarts;
     5601                           delete [] addRow;
     5602                           delete [] addElement;
     5603                           delete [] addCost;
     5604                         }
     5605#endif
    51385606                         // Set rowCounts to be diagonal block index for cleanup
    51395607                         rowCounts[kBlock] = jBlock;
     
    51415609                         // master
    51425610                         masterBlock = jBlock;
    5143                          master.loadProblem(*matrix, columnLower, columnUpper,
     5611                         masterModel.loadProblem(*matrix, columnLower, columnUpper,
    51445612                                            objective, rowLower, rowUpper);
    51455613                         if (optimizationDirection_ < 0.0) {
    5146                               double * obj = master.objective();
    5147                               int n = master.numberColumns();
     5614                              double * obj = masterModel.objective();
     5615                              int n = masterModel.numberColumns();
    51485616                              for (int i = 0; i < n; i++)
    51495617                                   obj[i] = - obj[i];
     
    51575625     delete [] whichBlock;
    51585626     delete [] blockInfo;
    5159      // For now master must have been defined (does not have to have rows)
    5160      assert (master.numberColumns());
    51615627     assert (masterBlock >= 0);
    5162      int numberMasterColumns = master.numberColumns();
     5628     int numberMasterColumns = masterModel.numberColumns();
     5629     masterModel.setStrParam(ClpProbName,"Master");
    51635630     // Overkill in terms of space
    51645631     int spaceNeeded = CoinMax(numberBlocks * (numberMasterColumns + 1),
     
    51675634     double * elementAdd = new double[spaceNeeded];
    51685635     spaceNeeded = numberBlocks;
    5169      int * rowAdd = new int[spaceNeeded+1];
     5636     int * rowAdd = new int[2*spaceNeeded+1]; // temp for block info
     5637     int * blockPrint = rowAdd+spaceNeeded+1;
    51705638     double * objective = new double[spaceNeeded];
    5171      int maxPass = 500;
     5639     int logLevel=handler_->logLevel();
     5640     //#define TEST_MODEL
     5641#ifdef TEST_MODEL
     5642     double goodValue=COIN_DBL_MAX;
     5643     ClpSimplex goodModel;
     5644     if (logLevel>3) {
     5645       // temp - create copy with master at front
     5646       const int * columnBack = model->coinBlock(masterBlock)->originalColumns();
     5647       int numberColumns2 = model->coinBlock(masterBlock)->numberColumns();
     5648       const int * rowBack = model->coinBlock(masterBlock)->originalRows();
     5649       int numberRows2 = model->coinBlock(masterBlock)->numberRows();
     5650       int * whichColumn = new int [numberColumns_+numberBlocks];
     5651       int * whichRow = new int [numberRows_];
     5652       int nColumn=0;
     5653       for (int iColumn = 0; iColumn < numberColumns2; iColumn++) {
     5654         int kColumn = columnBack[iColumn];
     5655         whichColumn[nColumn++]=kColumn;
     5656       }
     5657       int nRow=0;
     5658       for (int iRow = 0; iRow < numberRows2; iRow++) {
     5659         int kRow = rowBack[iRow];
     5660         whichRow[nRow++]=kRow;
     5661       }
     5662       for (iBlock = 0; iBlock < numberBlocks; iBlock++) {
     5663         int kBlock = rowCounts[iBlock];
     5664         const int * columnBack = model->coinBlock(kBlock)->originalColumns();
     5665         int numberColumns2 = model->coinBlock(kBlock)->numberColumns();
     5666         const int * rowBack = model->coinBlock(kBlock)->originalRows();
     5667         int numberRows2 = model->coinBlock(kBlock)->numberRows();
     5668         for (int iColumn = 0; iColumn < numberColumns2; iColumn++) {
     5669           int kColumn = columnBack[iColumn];
     5670           whichColumn[nColumn++]=kColumn;
     5671         }
     5672         for (int iRow = 0; iRow < numberRows2; iRow++) {
     5673           int kRow = rowBack[iRow];
     5674           whichRow[nRow++]=kRow;
     5675         }
     5676       }
     5677       ClpSimplex temp(this,nRow,whichRow,nColumn,whichColumn);
     5678       temp.writeMps("ordered.mps");
     5679       for (int i=numberMasterColumns;i<numberColumns_;i++)
     5680         whichColumn[i+numberBlocks]=i;
     5681       for (int i=0;i<numberMasterColumns;i++)
     5682         whichColumn[i]=i;
     5683       for (int i=0;i<numberBlocks;i++)
     5684         whichColumn[i+numberMasterColumns]=i+numberColumns_;
     5685       double * lower = new double [numberBlocks];
     5686       // Add extra variables
     5687       columnAdd[0] = 0;
     5688       for (int iBlock = 0; iBlock < numberBlocks; iBlock++) {
     5689         objective[iBlock] = 0.0;
     5690         lower[iBlock]=-COIN_DBL_MAX;
     5691         columnAdd[iBlock+1] = 0;
     5692       }
     5693       temp.addColumns(numberBlocks, lower, NULL, objective,
     5694                       columnAdd, rowAdd, elementAdd);
     5695       delete [] lower;
     5696       ClpSimplex temp2(&temp,nRow,whichRow,nColumn+numberBlocks,whichColumn);
     5697       goodModel=temp2;
     5698       goodModel.dual();
     5699       goodValue=goodModel.objectiveValue();
     5700#if 0
     5701       double * obj = goodModel.objective();
     5702       for (int i=numberMasterColumns;i<numberMasterColumns+numberBlocks;i++)
     5703         obj[i]=0.5;;
     5704       for (int i=numberMasterColumns+numberBlocks;i<goodModel.numberColumns();i++)
     5705         obj[i]*=0.5;
     5706       // fix solution
     5707       {
     5708         double * lower = goodModel.columnLower();
     5709         double * upper = goodModel.columnUpper();
     5710         double * solution = goodModel.primalColumnSolution();
     5711         for (int i=0;i<numberMasterColumns;i++) {
     5712           lower[i]=solution[i];
     5713           upper[i]=solution[i];
     5714         }
     5715         goodModel.scaling(0);
     5716         goodModel.dual();
     5717       }
     5718#endif
     5719       delete [] whichColumn;
     5720       delete [] whichRow;
     5721     }
     5722#endif
     5723     int maxPass = options.independentOption(2);
     5724     if (maxPass<2)
     5725       maxPass=100;
    51725726     int iPass;
    51735727     double lastObjective = -1.0e31;
    51745728     // Create columns for proposals
    5175      int numberMasterRows = master.numberRows();
    5176      master.resize(numberMasterColumns + numberBlocks, numberMasterRows);
     5729     int numberMasterRows = masterModel.numberRows();
     5730     //masterModel.resize(numberMasterColumns + numberBlocks, numberMasterRows);
    51775731     if (this->factorizationFrequency() == 200) {
    51785732          // User did not touch preset
    5179           master.defaultFactorizationFrequency();
     5733          masterModel.defaultFactorizationFrequency();
    51805734     } else {
    51815735          // make sure model has correct value
    5182           master.setFactorizationFrequency(this->factorizationFrequency());
    5183      }
    5184      master.setPerturbation(50);
     5736          masterModel.setFactorizationFrequency(this->factorizationFrequency());
     5737     }
     5738     masterModel.setPerturbation(50);
     5739     // temp bounds
     5740     if (0) {
     5741       printf("temp bounds\n");
     5742       double * lower = masterModel.columnLower();
     5743       double * upper = masterModel.columnUpper();
     5744       for (int i=0;i<numberMasterColumns;i++) {
     5745         lower[i]=CoinMax(lower[i],-1.0e8);
     5746         upper[i]=CoinMin(upper[i],1.0e8);
     5747       }
     5748     }
     5749     //printf("take out debound\n");
     5750     //master=*deBound(&master);
    51855751     // Arrays to say which block and when created
    51865752     int maximumRows = 2 * numberMasterColumns + 10 * numberBlocks;
    51875753     whichBlock = new int[maximumRows];
    51885754     int * when = new int[maximumRows];
     5755     // state of each problem (0 first or infeas, 1 feasible, add 2 if extra variables freed)
     5756     int * problemState = new int [numberBlocks];
     5757     memset(problemState,0,numberBlocks*sizeof(int));
    51895758     int numberRowsGenerated = numberBlocks;
     5759     // space for rhs modifications
     5760     double ** modification = new double * [numberBlocks];
    51905761     // Add extra variables
    5191      {
    5192           int iBlock;
    5193           columnAdd[0] = 0;
    5194           for (iBlock = 0; iBlock < numberBlocks; iBlock++) {
    5195                objective[iBlock] = 1.0;
    5196                columnAdd[iBlock+1] = 0;
    5197                when[iBlock] = -1;
    5198                whichBlock[iBlock] = iBlock;
    5199           }
    5200           master.addColumns(numberBlocks, NULL, NULL, objective,
    5201                             columnAdd, rowAdd, elementAdd);
    5202      }
    5203      std::cout << "Time to decompose " << CoinCpuTime() - time1 << " seconds" << std::endl;
     5762     columnAdd[0] = 0;
     5763     for (int iBlock = 0; iBlock < numberBlocks; iBlock++) {
     5764       int numberRows2=sub[iBlock].numberRows();
     5765       int numberColumns2 = sub[iBlock].numberColumns();
     5766       int numberTotal2 = numberRows2+numberColumns2;
     5767       modification[iBlock]=new double [2*numberTotal2+numberRows2];
     5768       double * save = modification[iBlock];
     5769       memcpy(save,sub[iBlock].rowLower(),numberRows2*sizeof(double));
     5770       save += numberRows2;
     5771       memcpy(save,sub[iBlock].columnLower(),numberColumns2*sizeof(double));
     5772       save += numberColumns2;
     5773       memcpy(save,sub[iBlock].rowUpper(),numberRows2*sizeof(double));
     5774       save += numberRows2;
     5775       memcpy(save,sub[iBlock].columnUpper(),numberColumns2*sizeof(double));
     5776       objective[iBlock] = 1.0;
     5777       columnAdd[iBlock+1] = 0;
     5778       when[iBlock] = -1;
     5779       whichBlock[iBlock] = iBlock;
     5780       if (logLevel<2) {
     5781         // modify printing
     5782         sub[iBlock].messagesPointer()->setDetailMessage(2,6);
     5783         sub[iBlock].messagesPointer()->setDetailMessage(2,14);
     5784         sub[iBlock].messagesPointer()->setDetailMessage(2,0);
     5785       }
     5786       // scaling
     5787       sub[iBlock].scaling(scalingFlag_);
     5788       //sub[iBlock].scaling(0);
     5789     }
     5790     masterModel.addColumns(numberBlocks, NULL, NULL, objective,
     5791                       columnAdd, rowAdd, elementAdd);
     5792     sprintf(generalPrint,"Time to decompose %.2f seconds",CoinCpuTime() - time1);
     5793     handler_->message(CLP_GENERAL, messages_)
     5794       << generalPrint
     5795       << CoinMessageEol;
     5796     int ixxxxxx=0;
     5797     if (intParam_[0]>=100000&&intParam_[0]<100999) {
     5798       ixxxxxx=intParam_[0]-100000;
     5799       printf("ixxxxxx %d\n",ixxxxxx);
     5800     }
     5801#ifdef ADD_ARTIFICIALS
     5802     double ** saveObjective = new double * [numberBlocks];
     5803     for (int i=0;i<numberBlocks;i++)
     5804       saveObjective[i]=CoinCopyOfArray(sub[i].objective(),
     5805                                        originalSubColumns[i]);
     5806     int iFudge=1;
     5807     int lowFudge=0;
     5808     int highFudge=0;
     5809#endif
     5810#define UNBOUNDED
     5811     if (ixxxxxx>0) {
     5812       for (iBlock = 0; iBlock < CoinMin(numberBlocks,ixxxxxx); iBlock++) {
     5813         ClpSimplex * temp = deBound(sub+iBlock);
     5814         sub[iBlock]=*temp;
     5815         delete temp;
     5816         delete [] modification[iBlock];
     5817         int numberRows2 = sub[iBlock].numberRows();
     5818         int numberColumns2 = sub[iBlock].numberColumns();
     5819         int numberTotal2 = numberRows2+numberColumns2;
     5820         modification[iBlock]=new double [2*numberTotal2+numberRows2];
     5821         double * save = modification[iBlock];
     5822         memcpy(save,sub[iBlock].rowLower(),numberRows2*sizeof(double));
     5823         save += numberRows2;
     5824         memcpy(save,sub[iBlock].columnLower(),numberColumns2*sizeof(double));
     5825         save += numberColumns2;
     5826         memcpy(save,sub[iBlock].rowUpper(),numberRows2*sizeof(double));
     5827         save += numberRows2;
     5828         memcpy(save,sub[iBlock].columnUpper(),numberColumns2*sizeof(double));
     5829       }
     5830     }
     5831#ifdef ABC_INHERIT
     5832     //AbcSimplex abcMaster;
     5833     //if (!this->abcState())
     5834     //setAbcState(1);
     5835     int numberCpu=CoinMin((this->abcState()&15),4);
     5836     AbcPthreadStuff threadInfo(numberCpu);
     5837     masterModel.setAbcState(this->abcState());
     5838     //AbcSimplex * tempMaster=masterModel.dealWithAbc(2,10,true);
     5839     //abcMaster=*tempMaster;
     5840     //delete tempMaster;
     5841     //abcMaster.startThreads(numberCpu);
     5842     //#define masterModel abcMaster
     5843#endif
     5844     double treatSubAsFeasible=1.0e-6;
     5845     int numberSubInfeasible=0;
     5846     bool canSkipSubSolve=false;
     5847     int numberProposals=999;
    52045848     for (iPass = 0; iPass < maxPass; iPass++) {
    5205           printf("Start of pass %d\n", iPass);
     5849          sprintf(generalPrint,"Start of pass %d", iPass);
     5850          handler_->message(CLP_GENERAL, messages_)
     5851            << generalPrint
     5852            << CoinMessageEol;
    52065853          // Solve master - may be unbounded
    5207           //master.scaling(0);
    5208           if (1) {
    5209                master.writeMps("yy.mps");
    5210           }
    5211           master.dual();
    5212           int problemStatus = master.status(); // do here as can change (delcols)
    5213           if (master.numberIterations() == 0 && iPass)
    5214                break; // finished
    5215           if (master.objectiveValue() < lastObjective + 1.0e-7 && iPass > 555)
    5216                break; // finished
    5217           lastObjective = master.objectiveValue();
     5854          //masterModel.scaling(0);
     5855          // get obj for debug
     5856          double objSum=masterModel.objectiveValue();
     5857          for (int i=0;i<numberBlocks;i++)
     5858            objSum += sub[i].objectiveValue();
     5859          //printf("objsum %g\n",objSum);
     5860          if (0) {
     5861               masterModel.writeMps("yy.mps");
     5862               masterModel.writeBasis("yy.bas",true,2);
     5863          }
     5864          {
     5865            // free up extra variables
     5866            double * lower = masterModel.columnLower();
     5867            int numberFreed=0;
     5868            for (int i=0;i<numberBlocks;i++) {
     5869              if (problemState[i]==1
     5870                  ) {
     5871                // ? need trust region ?
     5872                lower[i+numberMasterColumns]=-COIN_DBL_MAX;
     5873#if 0 //1 //ndef UNBOUNDED
     5874                lower[i+numberMasterColumns]=-1.0e10;
     5875#endif
     5876                //if (problemState[i]!=2) {
     5877                  numberFreed++;
     5878                  problemState[i]=3;
     5879                  //}
     5880              }
     5881            }
     5882            if (numberFreed)
     5883              lastObjective = -1.0e31;
     5884          }
     5885#ifdef TRY_NO_SCALING
     5886          masterModel.scaling(0);
     5887#endif
     5888#ifdef ABC_INHERIT
     5889          masterModel.dealWithAbc(0,0,true);
     5890#else
     5891          masterModel.dual();
     5892#endif
     5893          if ((maxPass==5000&&scalingFlag_)||(maxPass==4000&&!scalingFlag_)) {
     5894            int n=masterModel.numberIterations();
     5895            masterModel.scaling(0);
     5896            masterModel.primal();
     5897            masterModel.setNumberIterations(n+masterModel.numberIterations());
     5898            masterModel.scaling(scalingFlag_);
     5899          }
     5900          int masterStatus = masterModel.status(); // do here as can change
     5901          sprintf(generalPrint,"Pass %d objective %g change %g",
     5902                 iPass,masterModel.objectiveValue(),
     5903                 masterModel.objectiveValue()-lastObjective);
     5904          handler_->message(CLP_GENERAL, messages_)
     5905            << generalPrint
     5906            << CoinMessageEol;
     5907#ifndef UNBOUNDED
     5908          if(masterStatus==2) {
     5909            // unbounded
     5910            masterModel.writeMps("unbounded.mps");
     5911            // get primal feasible
     5912#ifdef ABC_INHERIT
     5913            masterModel.dealWithAbc(1,1,true);
     5914#else
     5915            masterModel.primal();
     5916#endif
     5917            const double * fullLower = columnLower();
     5918            const double * fullUpper = columnUpper();
     5919            double * lower = masterModel.columnLower();
     5920            double * upper = masterModel.columnUpper();
     5921            double * solution = masterModel.primalColumnSolution();
     5922            const int * columnBack =
     5923              model->coinBlock(masterBlock)->originalColumns();
     5924            if (!numberTimesUnbounded) {
     5925              for (int iColumn = 0; iColumn < numberMasterColumns; iColumn++) {
     5926                int kColumn = columnBack[iColumn];
     5927                double value = solution[iColumn];
     5928                double lowerValue=CoinMax(fullLower[kColumn],
     5929                                          CoinMin(value,fullUpper[kColumn])-trust);
     5930                lower[iColumn]=lowerValue;
     5931                double upperValue=CoinMin(fullUpper[kColumn],
     5932                                          CoinMax(value,fullLower[kColumn])+trust);
     5933                upper[iColumn]=upperValue;
     5934              }
     5935#ifdef TEST_MODEL
     5936              if (logLevel>3) {
     5937                //use ones from solved
     5938                const double * solutionGood = goodModel.primalColumnSolution();
     5939                for (int iColumn = 0; iColumn < numberMasterColumns; iColumn++) {
     5940                  double value = solutionGood[iColumn];
     5941                  lower[iColumn]=CoinMin(value,-trust);
     5942                  upper[iColumn]=CoinMax(value,trust);
     5943                }
     5944              }
     5945#endif
     5946            } else {
     5947              abort(); // probably can happen
     5948            }
     5949            numberTimesUnbounded++;
     5950#ifdef ABC_INHERIT
     5951            masterModel.dealWithAbc(0,0,true);
     5952#else
     5953            masterModel.dual();
     5954#endif
     5955            masterStatus = masterModel.status();
     5956            assert (!masterStatus);
     5957            masterModel.setNumberIterations(1); // so will continue
     5958          }
     5959          if (numberTimesUnbounded>1&&trust>0.0) {
     5960            assert (!masterStatus);
     5961            const double * fullLower = columnLower();
     5962            const double * fullUpper = columnUpper();
     5963            double * lower = masterModel.columnLower();
     5964            double * upper = masterModel.columnUpper();
     5965            //double * solution = masterModel.primalColumnSolution();
     5966            const int * columnBack =
     5967              model->coinBlock(masterBlock)->originalColumns();
     5968            int nTrusted=0;
     5969            for (int iColumn = 0; iColumn < numberMasterColumns; iColumn++) {
     5970              int kColumn = columnBack[iColumn];
     5971              if (lower[iColumn]>fullLower[kColumn]) {
     5972                if(masterModel.getColumnStatus(iColumn)!=atLowerBound) {
     5973                  lower[iColumn]=fullLower[kColumn];
     5974                } else {
     5975                  nTrusted++;
     5976                }
     5977              }
     5978              if (upper[iColumn]<fullUpper[kColumn]) {
     5979                if(masterModel.getColumnStatus(iColumn)!=atUpperBound) {
     5980                  upper[iColumn]=fullUpper[kColumn];
     5981                } else {
     5982                  nTrusted++;
     5983                }
     5984              }
     5985            }
     5986            if (nTrusted) {
     5987              sprintf(generalPrint,"%d at artificial bound",nTrusted);
     5988            } else {
     5989              sprintf(generalPrint,"All at natural bounds");
     5990              trust=0.0;
     5991            }
     5992            handler_->message(CLP_GENERAL2, messages_)
     5993            << generalPrint
     5994            << CoinMessageEol;
     5995          }
     5996#endif
     5997          if (!masterStatus) {
     5998            if (masterModel.numberIterations() == 0 && iPass ) {
     5999              if ((!numberSubInfeasible&&!numberProposals)||treatSubAsFeasible>1.0e-2||iPass>5555)
     6000                break; // finished
     6001              if (!numberProposals&&numberSubInfeasible) {
     6002                treatSubAsFeasible *= 2.0;
     6003                printf("Doubling sub primal tolerance to %g\n",treatSubAsFeasible);
     6004              } else {
     6005                treatSubAsFeasible *= 1.2;
     6006                printf("Increasing sub primal tolerance to %g\n",treatSubAsFeasible);
     6007              }
     6008              canSkipSubSolve=false;
     6009            } else if (!numberSubInfeasible) {
     6010              if (treatSubAsFeasible>1.0e-6) {
     6011                treatSubAsFeasible = CoinMax(0.9*treatSubAsFeasible,1.0e-6);
     6012                printf("Reducing sub primal tolerance to %g\n",treatSubAsFeasible);
     6013              }
     6014            }         
     6015            if (masterModel.objectiveValue() < lastObjective + 1.0e-7 && iPass > 5555)
     6016              break; // finished
     6017            lastObjective = masterModel.objectiveValue();
     6018          }
    52186019          // mark non-basic rows and delete if necessary
    52196020          int iRow;
    5220           numberRowsGenerated = master.numberRows() - numberMasterRows;
     6021          numberRowsGenerated = masterModel.numberRows() - numberMasterRows;
    52216022          for (iRow = 0; iRow < numberRowsGenerated; iRow++) {
    5222                if (master.getStatus(iRow + numberMasterRows) != ClpSimplex::basic)
     6023               if (masterModel.getStatus(iRow + numberMasterRows) != ClpSimplex::basic)
    52236024                    when[iRow] = iPass;
    52246025          }
    5225           if (numberRowsGenerated > maximumRows) {
     6026          if (numberRowsGenerated > maximumRows-numberBlocks) {
    52266027               // delete
    52276028               int numberKeep = 0;
     
    52296030               int * whichDelete = new int[numberRowsGenerated];
    52306031               for (iRow = 0; iRow < numberRowsGenerated; iRow++) {
    5231                     if (when[iRow] > iPass - 7) {
    5232                          // keep
    5233                          when[numberKeep] = when[iRow];
    5234                          whichBlock[numberKeep++] = whichBlock[iRow];
    5235                     } else {
    5236                          // delete
    5237                          whichDelete[numberDelete++] = iRow + numberMasterRows;
    5238                     }
    5239                }
     6032                 if (masterModel.getRowStatus(iRow+numberMasterRows)!= basic) {
     6033                   // keep
     6034                   when[numberKeep] = when[iRow];
     6035                   whichBlock[numberKeep++] = whichBlock[iRow];
     6036                 } else {
     6037                   // delete
     6038                   whichDelete[numberDelete++] = iRow + numberMasterRows;
     6039                 }
     6040               }
     6041               if (numberRowsGenerated-numberDelete > maximumRows-numberBlocks) {
     6042                 for (iRow = 0; iRow < numberRowsGenerated; iRow++) {
     6043                   if (when[iRow] > iPass - 7 ) {
     6044                     // keep
     6045                     when[numberKeep] = when[iRow];
     6046                     whichBlock[numberKeep++] = whichBlock[iRow];
     6047                   } else {
     6048                     // delete
     6049                     whichDelete[numberDelete++] = iRow + numberMasterRows;
     6050                   }
     6051                 }
     6052               }
    52406053               numberRowsGenerated -= numberDelete;
    5241                master.deleteRows(numberDelete, whichDelete);
     6054               masterModel.deleteRows(numberDelete, whichDelete);
    52426055               delete [] whichDelete;
    52436056          }
    5244           const double * primal = NULL;
     6057          double * primal = NULL;
    52456058          bool deletePrimal = false;
    5246           if (problemStatus == 0) {
    5247                primal = master.primalColumnSolution();
    5248           } else if (problemStatus == 2) {
    5249                // could do composite objective
    5250                primal = master.unboundedRay();
    5251                deletePrimal = true;
    5252                printf("The sum of infeasibilities is %g\n",
    5253                       master.sumPrimalInfeasibilities());
    5254           } else if (!master.numberRows()) {
     6059          if (masterStatus == 0) {
     6060               primal = masterModel.primalColumnSolution();
     6061          } else if (masterStatus == 2 && masterModel.numberRows()) {
     6062               // scale back ray (1.0e20?)
     6063               primal = masterModel.ray();
     6064               //deletePrimal = true;
     6065               sprintf(generalPrint,"The sum of dual infeasibilities is %g",
     6066                      masterModel.sumDualInfeasibilities());
     6067               handler_->message(CLP_GENERAL, messages_)
     6068                 << generalPrint
     6069                 << CoinMessageEol;
     6070#if 0
     6071               const double * primal2 = masterModel.primalColumnSolution();
     6072#ifndef UNBOUNDED
     6073               for (int i=0;i<numberMasterColumns;i++) {
     6074                 primal[i] = 1.0e-10*primal[i]+primal2[i];
     6075               }
     6076#else
     6077               double scaleFactor = innerProduct(primal,numberMasterColumns,primal);
     6078               double scaleFactor2 = innerProduct(primal+numberMasterColumns,
     6079                                                  numberBlocks,primal+numberMasterColumns);
     6080               if (scaleFactor&&!scaleFactor2) {
     6081                 scaleFactor = 1.0/sqrt(scaleFactor);
     6082                 for (int i=0;i<numberMasterColumns;i++) {
     6083                   primal[i] *= scaleFactor;
     6084                 }
     6085               } else {
     6086                 // treat as feasible
     6087                 if (scaleFactor)
     6088                   scaleFactor = 1.0e10/sqrt(scaleFactor);
     6089                 for (int i=0;i<numberMasterColumns;i++) {
     6090                   primal[i] = primal[i]*scaleFactor+primal2[i];
     6091                 }
     6092                 masterStatus=0;
     6093               }
     6094#endif
     6095#endif
     6096          } else if (!masterModel.numberRows()) {
    52556097               assert(!iPass);
    5256                primal = master.primalColumnSolution();
    5257                memset(master.primalColumnSolution(),
     6098               primal = masterModel.primalColumnSolution();
     6099               memset(masterModel.primalColumnSolution(),
    52586100                      0, numberMasterColumns * sizeof(double));
    52596101          } else {
    5260                abort();
    5261           }
     6102            printf("Master infeasible - sum %g\n",
     6103                   masterModel.sumPrimalInfeasibilities());
     6104            masterModel.setProblemStatus(0);
     6105            primal = masterModel.primalColumnSolution();
     6106            masterModel.writeMps("inf.mps");
     6107            //abort();
     6108          }
     6109#ifndef UNBOUNDED
     6110          if (masterStatus==2) {
     6111            // adjust variables with no elements
     6112            const int * columnLength = masterModel.matrix()->getVectorLengths();
     6113            const double * lower = masterModel.columnLower();
     6114            const double * upper = masterModel.columnUpper();
     6115            const double * obj = masterModel.objective();
     6116            for (int i=0;i<numberMasterColumns;i++) {
     6117              double value=primal[i];
     6118              if (!columnLength[i]) {
     6119                if (obj[i]<0.0)
     6120                  value += 1.0e10;
     6121                else if (obj[i]>0.0)
     6122                  value -= 1.0e10;
     6123              }
     6124              // make sure feasible
     6125              primal[i]=CoinMax(-1.0e10,CoinMin(1.0e10,value));
     6126              primal[i]=CoinMax(lower[i],CoinMin(upper[i],primal[i]));
     6127            }
     6128          }
     6129#endif
    52626130          // Create rhs for sub problems and solve
     6131          for (iBlock = 0; iBlock < numberBlocks; iBlock++) {
     6132#ifdef ADD_ARTIFICIALS
     6133            {
     6134              double * columnUpper2=sub[iBlock].columnUpper();
     6135              double * obj = sub[iBlock].objective();
     6136              int start=originalSubColumns[iBlock];
     6137              int numberColumns2=sub[iBlock].numberColumns();
     6138              if (iFudge>=lowFudge&&iFudge<=abs(highFudge)) {
     6139                for (int i=start;i<numberColumns2;i++)
     6140                  columnUpper2[i]=COIN_DBL_MAX;
     6141                memset(obj,0,originalSubColumns[iBlock]*sizeof(double));
     6142              } else {
     6143                for (int i=start;i<numberColumns2;i++)
     6144                  columnUpper2[i]=0.0;
     6145                memcpy(obj,saveObjective[iBlock],originalSubColumns[iBlock]*sizeof(double));
     6146                if (highFudge<0) {
     6147                  sub[iBlock].allSlackBasis(true);
     6148                }
     6149              }
     6150            }
     6151#endif
     6152            int numberRows2 = sub[iBlock].numberRows();
     6153            int numberColumns2 = sub[iBlock].numberColumns();
     6154            double * saveLower = modification[iBlock];
     6155            double * lower2 = sub[iBlock].rowLower();
     6156            memcpy(lower2, saveLower, numberRows2 * sizeof(double));
     6157            double * saveColumnLower = saveLower+numberRows2;
     6158            double * columnLower2 = sub[iBlock].columnLower();
     6159            memcpy(columnLower2,saveColumnLower,numberColumns2*sizeof(double));
     6160            double * saveUpper = saveColumnLower+numberColumns2;
     6161            double * upper2 = sub[iBlock].rowUpper();
     6162            memcpy(upper2, saveUpper, numberRows2 * sizeof(double));
     6163            double * saveColumnUpper = saveUpper+numberRows2;
     6164            double * columnUpper2 = sub[iBlock].columnUpper();
     6165            memcpy(columnUpper2,saveColumnUpper,numberColumns2*sizeof(double));
     6166            double * lastMod = saveColumnUpper+numberColumns2;
     6167#ifdef UNBOUNDED
     6168            if (masterStatus==2) {
     6169              for (int i=0;i<numberRows2;i++) {
     6170                if (lower2[i]>-COIN_DBL_MAX)
     6171                  lower2[i]=0.0;
     6172                if (upper2[i]<COIN_DBL_MAX)
     6173                  upper2[i]=0.0;
     6174              }
     6175              for (int i=0;i<numberColumns2;i++) {
     6176                if (columnLower2[i]>-COIN_DBL_MAX)
     6177                  columnLower2[i]=0.0;
     6178                if (columnUpper2[i]<COIN_DBL_MAX)
     6179                  columnUpper2[i]=0.0;
     6180              }
     6181            }
     6182#endif
     6183            // new rhs
     6184            double * rhs = sub[iBlock].dualRowSolution();
     6185            CoinZeroN(rhs, numberRows2);
     6186            first[iBlock]->times(primal, rhs);
     6187            for (int i = 0; i < numberRows2; i++) {
     6188              double value = rhs[i];
     6189              if (lower2[i] > -1.0e30)
     6190                lower2[i] -= value;
     6191              if (upper2[i] < 1.0e30)
     6192                upper2[i] -= value;
     6193            }
     6194            bool canSkip=false;
     6195            if (canSkipSubSolve) {
     6196              canSkip=true;
     6197              const double * rowSolution = sub[iBlock].primalRowSolution();
     6198              for (int i = 0; i < numberRows2; i++) {
     6199                double value = lastMod[i]-rhs[i];
     6200                if (fabs(value)>primalTolerance_) {
     6201                  // see if we can adjust
     6202                  double rowValue=rowSolution[i];
     6203                  if (rowValue<lower2[i]-primalTolerance_
     6204                      ||rowValue>upper2[i]+primalTolerance_) {
     6205                    canSkip=false;
     6206                    break;
     6207                  } else if (sub[iBlock].getRowStatus(i)!=basic) {
     6208                    canSkip=false;
     6209                    break;
     6210                  }
     6211                }
     6212              }
     6213            }
     6214            if (!canSkip) {
     6215              memcpy(lastMod,rhs,numberRows2*sizeof(double));
     6216            } else {
     6217              // mark
     6218              sub[iBlock].setSecondaryStatus(99);
     6219            }
     6220          }
     6221          canSkipSubSolve=true;
     6222          if (!iPass) {
     6223            // do first and then copy status?
     6224#ifdef TRY_NO_SCALING
     6225            sub[0].scaling(0);
     6226#endif
     6227#ifdef ABC_INHERIT
     6228            sub[0].setAbcState(numberCpu);
     6229            //sub[0].dealWithAbc(0,0,true);
     6230            sub[0].dealWithAbc(1,1,true);
     6231            sub[0].setAbcState(0);
     6232#else
     6233            sub[0].dual();
     6234#endif
     6235            if ((maxPass==5000&&scalingFlag_)||(maxPass==4000&&!scalingFlag_)) {
     6236              int n=sub[0].numberIterations();
     6237              sub[0].scaling(0);
     6238              sub[0].primal();
     6239              sub[0].setNumberIterations(n+sub[0].numberIterations());
     6240              sub[0].scaling(scalingFlag_);
     6241            }
     6242            int numberIterations=sub[0].numberIterations();
     6243            if (sub[0].problemStatus()) {
     6244              sub[0].primal();
     6245#if 0
     6246              sub[0].writeMps("first.mps");
     6247              sub[0].writeBasis("first.bas",true);
     6248              sub[0].readBasis("first.bas");
     6249              sub[0].primal();
     6250#endif
     6251              numberIterations+=sub[0].numberIterations();
     6252            }
     6253            sprintf(generalPrint,"First block - initial solve - %d iterations, objective %g",
     6254                   numberIterations,sub[0].objectiveValue());
     6255            handler_->message(CLP_GENERAL, messages_)
     6256              << generalPrint
     6257              << CoinMessageEol;
     6258            // copy status if same size
     6259            int numberRows2 = sub[0].numberRows();
     6260            int numberColumns2 = sub[0].numberColumns();
     6261            for (int iBlock = 1; iBlock < numberBlocks; iBlock++) {
     6262              int numberRows2a = sub[iBlock].numberRows();
     6263              int numberColumns2a = sub[iBlock].numberColumns();
     6264              if (numberRows2==numberRows2a&&
     6265                  numberColumns2==numberColumns2a) {
     6266                memcpy(sub[iBlock].primalColumnSolution(),
     6267                       sub[0].primalColumnSolution(),
     6268                       numberColumns2*sizeof(double));
     6269                memcpy(sub[iBlock].statusArray(),sub[0].statusArray(),
     6270                       numberRows2+numberColumns2);
     6271              }
     6272            }
     6273            // mark
     6274            sub[0].setSecondaryStatus(99);
     6275          }
     6276#ifdef ABC_INHERIT
     6277          if (numberCpu<2) {
     6278#endif
     6279            numberSubInfeasible=0;
     6280            for (iBlock = 0; iBlock < numberBlocks; iBlock++) {
     6281#ifdef TRY_NO_SCALING
     6282              sub[iBlock].scaling(0);
     6283#endif
     6284              if (sub[iBlock].secondaryStatus()!=99||false) {
     6285                //int ix=sub[iBlock].secondaryStatus();
     6286                int lastStatus = sub[iBlock].problemStatus();
     6287                // was do dual unless unbounded
     6288                double saveTolerance=sub[iBlock].primalTolerance();
     6289                if (lastStatus==0||!iPass) {
     6290                //if (lastStatus<2||!iPass) {
     6291                  //sub[iBlock].dual();
     6292                  sub[iBlock].primal();
     6293                  if (!sub[iBlock].isProvenOptimal()&&
     6294                      sub[iBlock].sumPrimalInfeasibilities()<treatSubAsFeasible) {
     6295                    printf("Block %d was feasible now has small infeasibility %g\n",iBlock,
     6296                           sub[iBlock].sumPrimalInfeasibilities());
     6297                    sub[iBlock].setPrimalTolerance(CoinMin(treatSubAsFeasible,1.0e-4));
     6298                    sub[iBlock].setCurrentPrimalTolerance(CoinMin(treatSubAsFeasible,1.0e-4));
     6299                    sub[iBlock].primal();
     6300                    sub[iBlock].setProblemStatus(0);
     6301                    problemState[iBlock] |= 4; // force actions
     6302                  }
     6303                  if ((maxPass==5000&&scalingFlag_)||(maxPass==4000&&!scalingFlag_)) {
     6304                    int n=sub[iBlock].numberIterations();
     6305                    sub[iBlock].scaling(0);
     6306                    sub[iBlock].primal();
     6307                    sub[iBlock].setNumberIterations(n+sub[iBlock].numberIterations());
     6308                    sub[iBlock].scaling(scalingFlag_);
     6309                  }
     6310                } else if (lastStatus==1) {
     6311                  // zero out objective
     6312                  double saveScale = sub[iBlock].infeasibilityCost();
     6313                  ClpObjective * saveObjective = sub[iBlock].objectiveAsObject();
     6314                  int numberColumns = sub[iBlock].numberColumns();
     6315                  ClpLinearObjective fake(NULL,numberColumns);
     6316                  sub[iBlock].setObjectivePointer(&fake);
     6317                  int saveOptions = sub[iBlock].specialOptions();
     6318                  sub[iBlock].setSpecialOptions(saveOptions|8192);
     6319                  sub[iBlock].primal();
     6320                  if ((maxPass==5000&&scalingFlag_)||(maxPass==4000&&!scalingFlag_)) {
     6321                    int n=sub[iBlock].numberIterations();
     6322                    sub[iBlock].scaling(0);
     6323                    sub[iBlock].primal();
     6324                    sub[iBlock].setNumberIterations(n+sub[iBlock].numberIterations());
     6325                    sub[iBlock].scaling(scalingFlag_);
     6326                  }
     6327                  sub[iBlock].setObjectivePointer(saveObjective);
     6328                  sub[iBlock].setInfeasibilityCost(saveScale);
     6329                  if (!sub[iBlock].isProvenOptimal()&&
     6330                      sub[iBlock].sumPrimalInfeasibilities()<treatSubAsFeasible) {
     6331                    printf("Block %d was infeasible now has small infeasibility %g\n",iBlock,
     6332                           sub[iBlock].sumPrimalInfeasibilities());
     6333                    sub[iBlock].setProblemStatus(0);
     6334                    sub[iBlock].setPrimalTolerance(CoinMin(treatSubAsFeasible,1.0e-4));
     6335                    sub[iBlock].setCurrentPrimalTolerance(CoinMin(treatSubAsFeasible,1.0e-4));
     6336                  }
     6337                  if (sub[iBlock].isProvenOptimal()) {
     6338                    sub[iBlock].primal();
     6339                    if ((maxPass==5000&&scalingFlag_)||(maxPass==4000&&!scalingFlag_)) {
     6340                      int n=sub[iBlock].numberIterations();
     6341                      sub[iBlock].scaling(0);
     6342                      sub[iBlock].primal();
     6343                      sub[iBlock].setNumberIterations(n+sub[iBlock].numberIterations());
     6344                      sub[iBlock].scaling(scalingFlag_);
     6345                    }
     6346                    if (!sub[iBlock].isProvenOptimal()) {
     6347                      printf("Block %d infeasible on second go has small infeasibility %g\n",iBlock,
     6348                           sub[iBlock].sumPrimalInfeasibilities());
     6349                      sub[iBlock].setProblemStatus(0);
     6350                    }
     6351                    problemState[iBlock] |= 4; // force actions
     6352                  } else {
     6353                    printf("Block %d still infeasible - sum %g - %d iterations\n",iBlock,
     6354                           sub[iBlock].sumPrimalInfeasibilities(),
     6355                           sub[iBlock].numberIterations());
     6356                    numberSubInfeasible++;
     6357                    if (!sub[iBlock].ray()) {
     6358                      printf("Block %d has no ray!\n",iBlock);
     6359                      sub[iBlock].primal();
     6360                      assert (sub[iBlock].ray()); // otherwise declare optimal
     6361                    }
     6362                  }
     6363                  sub[iBlock].setSpecialOptions(saveOptions);
     6364                } else {
     6365                  sub[iBlock].primal();
     6366                }
     6367                sub[iBlock].setPrimalTolerance(saveTolerance);
     6368                sub[iBlock].setCurrentPrimalTolerance(saveTolerance);
     6369                if (!sub[iBlock].isProvenOptimal()&&!sub[iBlock].isProvenPrimalInfeasible()) {
     6370                  printf("!!!Block %d has bad status %d\n",iBlock,sub[iBlock].problemStatus());
     6371                  sub[iBlock].primal(); // last go
     6372                }
     6373                //#define WRITE_ALL
     6374#ifdef WRITE_ALL
     6375                char name[20];
     6376                sprintf(name,"pass_%d_block_%d.mps",iPass,iBlock);
     6377                sub[iBlock].writeMps(name);
     6378                sprintf(name,"pass_%d_block_%d.bas",iPass,iBlock);
     6379                sub[iBlock].writeBasis(name,true);
     6380                if (sub[iBlock].problemStatus()==1) {
     6381                  sub[iBlock].readBasis(name);
     6382                  sub[iBlock].primal();
     6383                }
     6384#endif
     6385                //assert (!sub[iBlock].numberIterations()||ix!=99);
     6386              }
     6387            }
     6388#ifdef ABC_INHERIT
     6389          } else {
     6390            int iBlock=0;
     6391            while (iBlock<numberBlocks) {
     6392              if (sub[iBlock].secondaryStatus()!=99) {
     6393                int iThread;
     6394                threadInfo.waitParallelTask(1,iThread,true);
     6395#ifdef DETAIL_THREAD
     6396                printf("Starting block %d on thread %d\n",
     6397                       iBlock,iThread);
     6398#endif
     6399                threadInfo.startParallelTask(1,iThread,sub+iBlock);
     6400              }
     6401              iBlock++;
     6402            }
     6403            threadInfo.waitAllTasks();
     6404          }
     6405#endif
     6406          for (iBlock = 0; iBlock < numberBlocks; iBlock++) {
     6407            if (!iPass)
     6408              problemState[iBlock] |= 4; // force actions
     6409            // if state changed then fake number of iterations
     6410            if ((problemState[iBlock]&1)==0) {
     6411              // was infeasible
     6412              if (sub[iBlock].isProvenOptimal()) {
     6413                // say feasible and changed
     6414                problemState[iBlock] |= 1+4;
     6415              }
     6416            } else {
     6417              // was feasible
     6418              if (sub[iBlock].isProvenPrimalInfeasible()) {
     6419                // say infeasible and changed
     6420                problemState[iBlock] &= ~1;
     6421                problemState[iBlock] |= 4;
     6422              }
     6423            }
     6424            if (sub[iBlock].secondaryStatus()!=99) {
     6425              if (logLevel>1) {
     6426                sprintf(generalPrint,"Block %d - %d iterations, objective %g",
     6427                       iBlock,sub[iBlock].numberIterations(),
     6428                       sub[iBlock].objectiveValue());
     6429                handler_->message(CLP_GENERAL2, messages_)
     6430                  << generalPrint
     6431                  << CoinMessageEol;
     6432              }
     6433              if (sub[iBlock].problemStatus()&&sub[iBlock].algorithm()<0&&false) {
     6434                int numberRows2=sub[iBlock].numberRows();
     6435                double * ray = sub[iBlock].infeasibilityRay();
     6436                double * saveRay=CoinCopyOfArray(ray,numberRows2);
     6437#if 0
     6438                int directionOut = sub[iBlock].directionOut();
     6439                sub[iBlock].allSlackBasis(true);
     6440                sub[iBlock].dual();
     6441                printf("old dir %d new %d\n",directionOut,sub[iBlock].directionOut());
     6442#else
     6443                double * obj = sub[iBlock].objective();
     6444                int numberColumns2=sub[iBlock].numberColumns();
     6445                double * saveObj = CoinCopyOfArray(obj,numberColumns2);
     6446                memset(obj,0,numberColumns2*sizeof(double));
     6447                sub[iBlock].allSlackBasis(true);
     6448                sub[iBlock].primal();
     6449                memcpy(obj,saveObj,numberColumns2*sizeof(double));
     6450                delete [] saveObj;
     6451#endif
     6452                ray = sub[iBlock].infeasibilityRay();
     6453                for (int i=0;i<numberRows2;i++) {
     6454                  if (fabs(ray[i]-saveRay[i])>1.0e-4+1.0e20)
     6455                    printf("** diffray block %d row %d first %g second %g\n",
     6456                           iBlock,i,saveRay[i],ray[i]);
     6457                }
     6458                delete [] saveRay;
     6459              }
     6460            }
     6461          }
     6462          if (!iPass)
     6463            sub[0].setSecondaryStatus(0);
     6464          if (logLevel>2) {
     6465            for (iBlock = 0; iBlock < numberBlocks; iBlock++) {
     6466              printf("block %d obj %g thetaC %g\n",iBlock,sub[iBlock].objectiveValue(),masterModel.primalColumnSolution()[numberMasterColumns+iBlock]);
     6467            }
     6468          }
    52636469          rowAdd[0] = 0;
    5264           int numberProposals = 0;
     6470          numberProposals = 0;
    52656471          for (iBlock = 0; iBlock < numberBlocks; iBlock++) {
    5266                int numberRows2 = sub[iBlock].numberRows();
    5267                double * saveLower = new double [numberRows2];
    5268                double * lower2 = sub[iBlock].rowLower();
    5269                double * saveUpper = new double [numberRows2];
    5270                double * upper2 = sub[iBlock].rowUpper();
    5271                // new rhs
    5272                CoinZeroN(saveUpper, numberRows2);
    5273                first[iBlock]->times(primal, saveUpper);
    5274                for (int i = 0; i < numberRows2; i++) {
    5275                     double value = saveUpper[i];
    5276                     saveLower[i] = lower2[i];
    5277                     saveUpper[i] = upper2[i];
    5278                     if (saveLower[i] > -1.0e30)
    5279                          lower2[i] -= value;
    5280                     if (saveUpper[i] < 1.0e30)
    5281                          upper2[i] -= value;
    5282                }
    5283                sub[iBlock].dual();
    5284                memcpy(lower2, saveLower, numberRows2 * sizeof(double));
    5285                memcpy(upper2, saveUpper, numberRows2 * sizeof(double));
    5286                // get proposal
    5287                if (sub[iBlock].numberIterations() || !iPass) {
    5288                     double objValue = 0.0;
    5289                     int start = rowAdd[numberProposals];
    5290                     // proposal
    5291                     if (sub[iBlock].isProvenOptimal()) {
    5292                          const double * solution = sub[iBlock].dualRowSolution();
    5293                          first[iBlock]->transposeTimes(solution, elementAdd + start);
    5294                          for (int i = 0; i < numberRows2; i++) {
    5295                               if (solution[i] < -dualTolerance_) {
    5296                                    // at upper
    5297                                    assert (saveUpper[i] < 1.0e30);
    5298                                    objValue += solution[i] * saveUpper[i];
    5299                               } else if (solution[i] > dualTolerance_) {
    5300                                    // at lower
    5301                                    assert (saveLower[i] > -1.0e30);
    5302                                    objValue += solution[i] * saveLower[i];
    5303                               }
    5304                          }
    5305 
    5306                          // See if cuts off and pack down
    5307                          int number = start;
    5308                          double infeas = objValue;
    5309                          double smallest = 1.0e100;
    5310                          double largest = 0.0;
    5311                          for (int i = 0; i < numberMasterColumns; i++) {
    5312                               double value = elementAdd[start+i];
    5313                               if (fabs(value) > 1.0e-15) {
    5314                                    infeas -= primal[i] * value;
    5315                                    smallest = CoinMin(smallest, fabs(value));
    5316                                    largest = CoinMax(largest, fabs(value));
    5317                                    columnAdd[number] = i;
    5318                                    elementAdd[number++] = -value;
    5319                               }
    5320                          }
    5321                          columnAdd[number] = numberMasterColumns + iBlock;
    5322                          elementAdd[number++] = -1.0;
    5323                          // if elements large then scale?
    5324                          //if (largest>1.0e8||smallest<1.0e-8)
    5325                          printf("For subproblem %d smallest - %g, largest %g - infeas %g\n",
    5326                                 iBlock, smallest, largest, infeas);
    5327                          if (infeas < -1.0e-6 || !iPass) {
    5328                               // take
    5329                               objective[numberProposals] = objValue;
    5330                               rowAdd[++numberProposals] = number;
    5331                               when[numberRowsGenerated] = iPass;
    5332                               whichBlock[numberRowsGenerated++] = iBlock;
    5333                          }
    5334                     } else if (sub[iBlock].isProvenPrimalInfeasible()) {
    5335                          // use ray
    5336                          const double * solution = sub[iBlock].infeasibilityRay();
    5337                          first[iBlock]->transposeTimes(solution, elementAdd + start);
    5338                          for (int i = 0; i < numberRows2; i++) {
    5339                               if (solution[i] < -dualTolerance_) {
    5340                                    // at upper
    5341                                    assert (saveUpper[i] < 1.0e30);
    5342                                    objValue += solution[i] * saveUpper[i];
    5343                               } else if (solution[i] > dualTolerance_) {
    5344                                    // at lower
    5345                                    assert (saveLower[i] > -1.0e30);
    5346                                    objValue += solution[i] * saveLower[i];
    5347                               }
    5348                          }
    5349                          // See if good infeas and pack down
    5350                          int number = start;
    5351                          double infeas = objValue;
    5352                          double smallest = 1.0e100;
    5353                          double largest = 0.0;
    5354                          for (int i = 0; i < numberMasterColumns; i++) {
    5355                               double value = elementAdd[start+i];
    5356                               if (fabs(value) > 1.0e-15) {
    5357                                    infeas -= primal[i] * value;
    5358                                    smallest = CoinMin(smallest, fabs(value));
    5359                                    largest = CoinMax(largest, fabs(value));
    5360                                    columnAdd[number] = i;
    5361                                    elementAdd[number++] = -value;
    5362                               }
    5363                          }
    5364                          // if elements large or small then scale?
    5365                          //if (largest>1.0e8||smallest<1.0e-8)
    5366                          printf("For subproblem ray %d smallest - %g, largest %g - infeas %g\n",
    5367                                 iBlock, smallest, largest, infeas);
    5368                          if (infeas < -1.0e-6) {
    5369                               // take
    5370                               objective[numberProposals] = objValue;
    5371                               rowAdd[++numberProposals] = number;
    5372                               when[numberRowsGenerated] = iPass;
    5373                               whichBlock[numberRowsGenerated++] = iBlock;
    5374                          }
    5375                     } else {
    5376                          abort();
    5377                     }
    5378                }
    5379                delete [] saveLower;
    5380                delete [] saveUpper;
    5381           }
    5382           if (deletePrimal)
    5383                delete [] primal;
     6472            int numberRows2 = sub[iBlock].numberRows();
     6473            int numberColumns2 = sub[iBlock].numberColumns();
     6474            double * saveLower = modification[iBlock];
     6475            double * lower2 = sub[iBlock].rowLower();
     6476            double * saveUpper = saveLower+numberRows2+numberColumns2;
     6477            double * upper2 = sub[iBlock].rowUpper();
     6478            int typeRun=sub[iBlock].secondaryStatus();
     6479            sub[iBlock].setSecondaryStatus(0);
     6480            if (typeRun!=99) {
     6481              if (0) {
     6482                double objValue = 0.0;
     6483                const double * solution = sub[iBlock].dualRowSolution();
     6484                for (int i = 0; i < numberRows2; i++) {
     6485                  if (solution[i] < -dualTolerance_) {
     6486                    // at upper
     6487                    assert (saveUpper[i] < 1.0e30);
     6488                    objValue += solution[i] * upper2[i];
     6489                  } else if (solution[i] > dualTolerance_) {
     6490                    // at lower
     6491                    assert (saveLower[i] > -1.0e30);
     6492                    objValue += solution[i] * lower2[i];
     6493                  }
     6494                }
     6495                //printf("obj %g\n",objValue);
     6496              }
     6497              // temp
     6498              if (sub[iBlock].isProvenPrimalInfeasible()&&!sub[iBlock].numberIterations())
     6499                problemState[iBlock] |= 4;
     6500              // get proposal
     6501              if (sub[iBlock].numberIterations()||(problemState[iBlock]&4)!=0) {
     6502                double objValue = 0.0;
     6503                int start = rowAdd[numberProposals];
     6504                // proposal
     6505                if (sub[iBlock].isProvenOptimal()) {
     6506                  double * solution = sub[iBlock].dualRowSolution();
     6507                  first[iBlock]->transposeTimes(solution, elementAdd + start);
     6508                  for (int i = 0; i < numberRows2; i++) {
     6509                    if (sub[iBlock].getRowStatus(i)==basic)
     6510                      solution[i]=0.0;
     6511                    if (saveUpper[i]>saveLower[i]) {
     6512                      if (sub[iBlock].getRowStatus(i)==atUpperBound)
     6513                        objValue += solution[i] * saveUpper[i];
     6514                      else
     6515                        objValue += solution[i] * saveLower[i];
     6516                    } else {
     6517                      // fixed
     6518                      objValue += solution[i] * saveLower[i];
     6519                    }
     6520                  }
     6521                  const double * dj = sub[iBlock].dualColumnSolution();
     6522                  const double * columnLower = sub[iBlock].columnLower();
     6523                  const double * columnUpper = sub[iBlock].columnUpper();
     6524                  double objValue2=0.0;
     6525                  int numberColumns2 = sub[iBlock].numberColumns();
     6526                  for (int i=0;i<numberColumns2;i++) {
     6527                    if (logLevel>2) {
     6528                      if (sub[iBlock].getColumnStatus(i)!=basic&&
     6529                          fabs(sub[iBlock].primalColumnSolution()[i])>1.0e-5)
     6530                        printf("zz %d has value %g\n",i,sub[iBlock].primalColumnSolution()[i]);
     6531                    }
     6532                    if (sub[iBlock].getColumnStatus(i)==isFixed) {
     6533                      objValue2+=columnLower[i]*dj[i];
     6534                    } else if (sub[iBlock].getColumnStatus(i)==atLowerBound) {
     6535                      objValue2+=columnLower[i]*dj[i];
     6536                    } else if (sub[iBlock].getColumnStatus(i)==atUpperBound) {
     6537                      objValue2+=columnUpper[i]*dj[i];
     6538                    }
     6539                  }
     6540#if 1
     6541                  double objValue3=0.0;
     6542                  const double * cost = sub[iBlock].objective();
     6543                  for (int i=0;i<numberColumns2;i++) {
     6544                    double value=dj[i]-cost[i];
     6545                    if (sub[iBlock].getColumnStatus(i)==isFixed) {
     6546                      objValue3+=columnLower[i]*value;
     6547                    } else if (sub[iBlock].getColumnStatus(i)==atLowerBound) {
     6548                      objValue3+=columnLower[i]*value;
     6549                    } else if (sub[iBlock].getColumnStatus(i)==atUpperBound) {
     6550                      objValue3+=columnUpper[i]*value;
     6551                    }
     6552                  }
     6553#endif
     6554                  // recompute
     6555                  if (logLevel>3) {
     6556                    printf("objValue %g from lp %g, obj2 %g, obj3 %g\n",
     6557                           objValue,sub[iBlock].objectiveValue(),
     6558                           objValue2,objValue3);
     6559                  }
     6560                  objValue += objValue2;
     6561                  //objValue=sub[iBlock].objectiveValue();
     6562                  // See if cuts off and pack down
     6563                  int number = start;
     6564                  double infeas = -objValue;
     6565                  double smallest = 1.0e100;
     6566                  double largest = 0.0;
     6567                  for (int i = 0; i < numberMasterColumns; i++) {
     6568                    double value = elementAdd[start+i];
     6569                    if (fabs(value) > 1.0e-12) {
     6570                      infeas += primal[i] * value;
     6571                      smallest = CoinMin(smallest, fabs(value));
     6572                      largest = CoinMax(largest, fabs(value));
     6573                      columnAdd[number] = i;
     6574                      elementAdd[number++] = -value;
     6575                    }
     6576                  }
     6577                 
     6578                  infeas += primal[numberMasterColumns+iBlock];
     6579                  columnAdd[number] = numberMasterColumns + iBlock;
     6580                  elementAdd[number++] = -1.0;
     6581                  // if elements large then scale?
     6582                  if (largest>1.0e8||smallest<1.0e-8) {
     6583                    sprintf(generalPrint,"For subproblem %d smallest - %g, largest %g - infeas %g",
     6584                           iBlock, smallest, largest, infeas);
     6585                    handler_->message(CLP_GENERAL2, messages_)
     6586                      << generalPrint
     6587                      << CoinMessageEol;
     6588                    if (smallest<1.0e-12*largest) {
     6589                      sprintf(generalPrint,"Removing small elements");
     6590                      handler_->message(CLP_GENERAL2, messages_)
     6591                        << generalPrint
     6592                        << CoinMessageEol;
     6593                      double target=1.0e-12*largest;
     6594                      smallest=largest;
     6595                      int number2=number-1;
     6596                      number = start;
     6597                      for (int i = start; i < number2; i++) {
     6598                        double value = elementAdd[i];
     6599                        if (fabs(value) > target) {
     6600                          smallest=CoinMin(smallest,fabs(value));
     6601                          columnAdd[number] = columnAdd[i];
     6602                          elementAdd[number++] = value;
     6603                        }
     6604                      }
     6605                      columnAdd[number] = numberMasterColumns + iBlock;
     6606                      elementAdd[number++] = -1.0;
     6607                    }
     6608                  }
     6609                  // if smallest >1.0 then scale
     6610                  if ((smallest>1.0e6||fabs(objValue)>0.01*largeValue_)
     6611                      &&number>start+1) {
     6612                    double scale = 1.0/smallest;
     6613                    if (fabs(scale*objValue)>0.01*largeValue_) {
     6614                      printf("** scale before obj scale %g\n",scale);
     6615                      scale =(0.01*largeValue_)/fabs(objValue);
     6616                    }
     6617                    printf("** scale %g infeas %g\n",scale,infeas);
     6618                    objValue *= scale;
     6619                    infeas = -objValue;
     6620                    for (int i = start; i < number-1; i++) {
     6621                      double value = elementAdd[i]*scale;
     6622                      elementAdd[i]=value;
     6623                      int iColumn=columnAdd[i];
     6624                      infeas -= primal[iColumn] * value;
     6625                    }
     6626                    elementAdd[number-1]*=scale;
     6627                    infeas += primal[numberMasterColumns+iBlock]*scale;
     6628                    printf("** new infeas %g - scales to %g\n",infeas,infeas/scale);
     6629                  }
     6630                  if (infeas < -1.0e-6 || (problemState[iBlock]&4)!=0) {
     6631                    // take
     6632                    // double check infeasibility
     6633                    if (logLevel>3)
     6634                      printf("objValue %g objectiveValue() %g\n",
     6635                             objValue,sub[iBlock].objectiveValue());
     6636                    double sum=0.0;
     6637                    for (int i=start;i<number;i++) {
     6638                      int iColumn=columnAdd[i];
     6639                      sum  += primal[iColumn]*elementAdd[i];
     6640                    }
     6641                    if (logLevel>3)
     6642                      printf("Sum %g rhs %g\n",sum,-objValue);
     6643                    if (logLevel>1)
     6644                      printf("Cut for block %d has %d elements\n",iBlock,number-1-start);
     6645                    blockPrint[numberProposals]=iBlock;
     6646                    objective[numberProposals] = -objValue;
     6647                    rowAdd[++numberProposals] = number;
     6648                    when[numberRowsGenerated] = iPass;
     6649                    whichBlock[numberRowsGenerated++] = iBlock;
     6650                  }
     6651                } else if (sub[iBlock].isProvenPrimalInfeasible()) {
     6652                  // use ray
     6653                  double * solution = sub[iBlock].infeasibilityRay();
     6654                  if (0) {
     6655                    double trueOffset=0.0;
     6656                    int numberRows=sub[iBlock].numberRows();
     6657                    int numberColumns=sub[iBlock].numberColumns();
     6658                    double * farkas = new double [CoinMax(2*numberColumns+numberRows,numberMasterColumns)];
     6659                    double * bound = farkas + numberColumns;
     6660                    double * effectiveRhs = bound + numberColumns;
     6661                    // get ray as user would
     6662                    double * ray = solution; //sub[iBlock].infeasibilityRay();
     6663                    // get farkas row
     6664                    memset(farkas,0,(2*numberColumns+numberRows)*sizeof(double));
     6665                    // Looks to me as if ray should be flipped according to mosek
     6666                    sub[iBlock].clpMatrix()->transposeTimes(-1.0,ray,farkas);
     6667                    // now farkas has A_T_y
     6668                    // Put nonzero bounds in bound
     6669                    const double * columnLower = sub[iBlock].columnLower();
     6670                    const double * columnUpper = sub[iBlock].columnUpper();
     6671                    int numberBad=0;
     6672                    // For sum in mosek
     6673                    double ySum=0.0;
     6674                    for (int i=0;i<numberColumns;i++) {
     6675                      double value = farkas[i];
     6676                      double boundValue=0.0;
     6677                      if (sub[iBlock].getStatus(i)==ClpSimplex::basic) {
     6678                        // treat as zero if small
     6679                        if (fabs(value)<1.0e-8) {
     6680                          value=0.0;
     6681                          farkas[i]=0.0;
     6682                        }
     6683                        if (value) {
     6684                          //printf("basic %d direction %d farkas %g\n",
     6685                          //       i,sub[iBlock].directionOut(),value);
     6686                          if (value<0.0)
     6687                            boundValue=columnLower[i];
     6688                          else
     6689                            boundValue=columnUpper[i];
     6690                        }
     6691                      } else if (fabs(value)>1.0e-8) {
     6692                        if (value<0.0)
     6693                          boundValue=columnLower[i];
     6694                        else
     6695                          boundValue=columnUpper[i];
     6696                        if (fabs(boundValue)>1.0e12&&fabs(value)<1.0e-8) {
     6697                          boundValue=0.0;
     6698                          value=0.0;
     6699                        }
     6700                      } else {
     6701                        value=0.0;
     6702                      }
     6703                      if (fabs(boundValue)>1.0e20) {
     6704                        numberBad++;
     6705                        boundValue=0.0;
     6706                        value=0.0;
     6707                        farkas[i]=0.0;
     6708                      }
     6709                      // mosek way
     6710                      // A_T_y + s_x_l -s_x_u == 0
     6711                      // So if value >0 s_x_l->0 s_x_u->value
     6712                      // otherwise s_x_l->-value, s_x_u->0
     6713                      double s_x_l=0.0;
     6714                      double s_x_u=0.0;
     6715                      if (value>0)
     6716                        s_x_u=value;
     6717                      else
     6718                        s_x_l=-value;
     6719                      ySum += columnLower[i]*s_x_l;
     6720                      ySum -= columnUpper[i]*s_x_u;
     6721                      bound[i]=boundValue;
     6722                    }
     6723                    const double * rowLower = sub[iBlock].rowLower();
     6724                    const double * rowUpper = sub[iBlock].rowUpper();
     6725                    //int pivotRow = sub[iBlock].spareIntArray_[3];
     6726                    //bool badPivot=pivotRow<0;
     6727                    for (int i=0;i<numberRows;i++) {
     6728                      double value = ray[i];
     6729                      double rhsValue=0.0;
     6730                      if (sub[iBlock].getRowStatus(i)==ClpSimplex::basic) {
     6731                        // treat as zero if small
     6732                        if (fabs(value)<1.0e-7) {
     6733                          value=0.0;
     6734                        }
     6735                        if (value) {
     6736                          //printf("row basic %d direction %d ray %g\n",
     6737                          //       i,sub[iBlock].directionOut(),value);
     6738                          if (value<0.0)
     6739                            rhsValue=rowLower[i];
     6740                          else
     6741                            rhsValue=rowUpper[i];
     6742                        }
     6743                      } else if (fabs(value)>1.0e-10) {
     6744                        if (value<0.0)
     6745                          rhsValue=rowLower[i];
     6746                        else
     6747                          rhsValue=rowUpper[i];
     6748                      } else {
     6749                        value=0.0;
     6750                      }
     6751                      if (fabs(rhsValue)>1.0e20) {
     6752                        numberBad++;
     6753                        value=0.0;
     6754                      }
     6755                      ray[i]=value;
     6756                      if (!value)
     6757                        rhsValue=0.0;
     6758                      // for mosek flip value back
     6759                      double yvalue = -value;
     6760                      // -y + s_c_l - s_c_u==0
     6761                      double s_c_l=0.0;
     6762                      double s_c_u=0.0;
     6763                      if (yvalue>0)
     6764                        s_c_l=yvalue;
     6765                      else
     6766                        s_c_u=-yvalue;
     6767                      ySum += rowLower[i]*s_c_l;
     6768                      ySum -= rowUpper[i]*s_c_u;
     6769                      effectiveRhs[i]=rhsValue;
     6770                      if (fabs(effectiveRhs[i])>1.0e10)
     6771                        printf("Large rhs row %d %g\n",
     6772                               i,effectiveRhs[i]);
     6773                    }
     6774                    sub[iBlock].clpMatrix()->times(-1.0,bound,effectiveRhs);
     6775                    double bSum=0.0;
     6776                    for (int i=0;i<numberRows;i++) {
     6777                      bSum += effectiveRhs[i]*ray[i];
     6778                      if (fabs(effectiveRhs[i])>1.0e10)
     6779                        printf("Large rhs row %d %g after\n",
     6780                               i,effectiveRhs[i]);
     6781                    }
     6782                    if (logLevel>1)
     6783                      printf("Block %d Mosek user manual wants %g to be positive so bSum should be negative %g\n",
     6784                             iBlock,ySum,bSum);
     6785                    if (numberBad||bSum>1.0e-6) {
     6786                      printf("Bad infeasibility ray %g  - %d bad\n",
     6787                             bSum,numberBad);
     6788                    } else {
     6789                      //printf("Good ray - infeasibility %g\n",
     6790                      //     -bSum);
     6791                    }
     6792                    /*
     6793                      wanted cut is
     6794                      plus or minus! (temp2 * x - temp2 *x_bar) <= bSum
     6795                      first[iBlock]->transposeTimes(ray, temp2);
     6796                     */
     6797                    memset(farkas,0,numberColumns*sizeof(double));
     6798                    first[iBlock]->transposeTimes(ray,farkas);
     6799                    double offset=0.0;
     6800                    const double * masterSolution = masterModel.primalColumnSolution();
     6801                    for (int i=0;i<numberMasterColumns;i++) {
     6802                      double value = farkas[i];
     6803                      if (fabs(value)>1.0e-9) {
     6804                        offset += value*masterSolution[i];
     6805                        if (logLevel>2)
     6806                          printf("(%d,%g) ",i,value);
     6807                      } else {
     6808                        farkas[i]=0.0;
     6809                      }
     6810                    }
     6811                    trueOffset = bSum+offset;
     6812                    if (sub[iBlock].algorithm()>0)
     6813                      trueOffset *= 1.0e-5;
     6814                    if (logLevel>2)
     6815                      printf(" - offset %g - ? rhs of %g\n",offset,trueOffset);
     6816                    //delete [] ray;
     6817                    delete [] farkas;
     6818                  }
     6819                  // if primal then scale
     6820                  if (sub[iBlock].algorithm()>0) {
     6821                    for (int i = 0; i < numberRows2; i++)
     6822                      solution[i] = -1.0e-5*solution[i];
     6823                  } else {
     6824                    for (int i = 0; i < numberRows2; i++)
     6825                      solution[i] = -solution[i];
     6826                  }
     6827                  first[iBlock]->transposeTimes(solution, elementAdd + start);
     6828                  for (int i = 0; i < numberRows2; i++)
     6829                    solution[i] = -solution[i];
     6830                  for (int i = 0; i < numberRows2; i++) {
     6831                    if (sub[iBlock].getRowStatus(i)==basic && fabs(solution[i])<1.0e-7)
     6832                      solution[i]=0.0;
     6833                    if (solution[i] > dualTolerance_) {
     6834                      // at upper
     6835                      if (saveUpper[i] > 1.0e20)
     6836                        solution[i]=0.0;
     6837                      objValue += solution[i] * saveUpper[i];
     6838                    } else if (solution[i] < -dualTolerance_) {
     6839                      // at lower
     6840                      if (saveLower[i] < -1.0e20);
     6841                        solution[i]=0.0;
     6842                      objValue += solution[i] * saveLower[i];
     6843                    } else {
     6844                      solution[i]=0.0;
     6845                    }
     6846                  }
     6847                  //objValue=-objValue;
     6848                  {
     6849                    int numberColumns2=sub[iBlock].numberColumns();
     6850                    double * temp = new double[numberColumns2];
     6851                    memset(temp,0,numberColumns2*sizeof(double));
     6852                    sub[iBlock].clpMatrix()->transposeTimes(-1.0,solution, temp);
     6853                    double loX=0.0;
     6854                    double upX=0.0;
     6855                    const double * lower = sub[iBlock].columnLower();
     6856                    const double * upper = sub[iBlock].columnUpper();
     6857                    const double * primal = sub[iBlock].primalColumnSolution();
     6858                    for (int i=0;i<numberColumns2;i++) {
     6859                      double value = temp[i];
     6860                      if (sub[iBlock].getColumnStatus(i)==basic&&
     6861                          fabs(value)<1.0e-7)
     6862                        value=0.0;
     6863                      if (logLevel>2) {
     6864                        //if (sub[iBlock].getColumnStatus(i)!=basic&&
     6865                        //  primal[i]>1.0e-5)
     6866                        //printf("zz_inf %d has value %g\n",i,primal[i]);
     6867                      }
     6868                      if (sub[iBlock].getStatus(i)==atLowerBound||
     6869                          sub[iBlock].getStatus(i)==isFixed) {
     6870                        loX += lower[i]*value;
     6871                      } else if (sub[iBlock].getStatus(i)==atUpperBound) {
     6872                        upX += upper[i]*value;
     6873                      } else if (sub[iBlock].getStatus(i)==basic) {
     6874                        double value2 = primal[i]*value;
     6875                        if (logLevel>2) {
     6876                          if (fabs(value2)>1.0e-3)
     6877                            printf("Basic %d arrayval %g primal %g bounds %g %g\n",
     6878                                   i,value,primal[i],
     6879                                   lower[i],upper[i]);
     6880                        }
     6881                        if (value<0.0) {
     6882                          assert (primal[i]<lower[i]);
     6883                          value2 = value*lower[i];
     6884                        } else if (value>0.0) {
     6885                          assert (primal[i]>upper[i]);
     6886                          if (primal[i]-upper[i]<1.0e-3) {
     6887                            if (logLevel>2)
     6888                              printf("small diff %g\n",primal[i]-upper[i]);
     6889                            //handler_->message(CLP_GENERAL2, messages_)
     6890                            //<< generalPrint
     6891                            //<< CoinMessageEol;
     6892                            //value=0.0;
     6893                            //elementAdd[start+i]=0.0;
     6894                          }
     6895                          value2 = value*upper[i];
     6896                        }
     6897                        loX += value2;
     6898                      }
     6899                    }
     6900                    objValue += loX+upX;
     6901                    if (logLevel>2)
     6902                      printf("Inf Offsets %g %g - new Objvalue %g\n",loX,upX,objValue);
     6903#define OBJ_OFFSET 0
     6904#if OBJ_OFFSET==1
     6905                    objValue -= loX+upX;
     6906                    objValue -= loX+upX;
     6907#elif OBJ_OFFSET==2
     6908                    objValue -= loX+upX;
     6909                    objValue = -objValue;
     6910                    objValue += loX+upX;
     6911#elif OBJ_OFFSET==3
     6912                    objValue -= loX+upX;
     6913                    objValue = -objValue;
     6914                    objValue -= loX+upX;
     6915#endif
     6916                    if (iBlock==-3) {
     6917                       ClpSimplex * temp = deBound(sub+iBlock);
     6918                       //temp->allSlackBasis();
     6919                       temp->primal();
     6920                       // use ray
     6921                       double * solution = temp->infeasibilityRay();
     6922                       int numberRows2=temp->numberRows();
     6923                       // bug somewhere - if primal then flip
     6924                       if (temp->algorithm_>0) {
     6925                         for (int i = 0; i < numberRows2; i++)
     6926                           solution[i] = - 1.0e-5 * solution[i];
     6927                       }
     6928                       double objValue7=0.0;
     6929                       const double * lower = temp->rowLower();
     6930                       const double * upper = temp->rowUpper();
     6931                       for (int i = 0; i < numberRows2; i++) {
     6932                         if (solution[i] < -dualTolerance_) {
     6933                           // at upper
     6934                           assert (upper[i] < 1.0e30);
     6935                           if (i<sub[iBlock].numberRows())
     6936                             objValue7 += solution[i] * saveUpper[i];
     6937                           else
     6938                             objValue7 += solution[i] * upper[i];
     6939                         } else if (solution[i] > dualTolerance_) {
     6940                           // at lower
     6941                           assert (lower[i] > -1.0e30);
     6942                           if (i<sub[iBlock].numberRows())
     6943                             objValue7 += solution[i] * saveLower[i];
     6944                           else
     6945                             objValue7 += solution[i] * lower[i];
     6946                         }
     6947                       }
     6948                       sprintf(generalPrint,"new objValue %g - old %g",objValue7,objValue);
     6949                       handler_->message(CLP_GENERAL2, messages_)
     6950                         << generalPrint
     6951                         << CoinMessageEol;
     6952                       //objValue=-objValue7;
     6953                       //loX=1.0e-2*objValue;
     6954                       //first[iBlock]->transposeTimes(solution, elementAdd + start);
     6955                       double * temp2 = new double [numberMasterColumns];
     6956                       memset(temp2,0,numberMasterColumns*sizeof(double));
     6957                       first[iBlock]->transposeTimes(solution, temp2);
     6958                       double * temp3 = elementAdd+start;
     6959                       for (int i=0;i<numberMasterColumns;i++) {
     6960                         if (fabs(temp2[i]-temp3[i])>1.0e-4)
     6961                           printf("** %d bound el %g nobound %g\n",i,temp3[i],temp2[i]);
     6962                       }
     6963                       memcpy(temp3,temp2,numberMasterColumns*sizeof(double));
     6964                       delete [] temp2;
     6965                       delete temp;
     6966                    }
     6967                    // relax slightly
     6968                    objValue += 1.0e-9*fabs(loX);
     6969                    delete [] temp;
     6970                  }
     6971                  delete [] solution;
     6972                  // See if good infeas and pack down (signs on infeas,value changed)
     6973                  int number = start;
     6974                  //printf("Not changing objValue from %g to %g\n",objValue,trueOffset);
     6975                  //printf("Changing objValue from %g to %g\n",objValue,trueOffset);
     6976                  //objValue=trueOffset;
     6977                  double infeas = objValue;
     6978                  double smallest = 1.0e100;
     6979                  double largest = 0.0;
     6980                  for (int i = 0; i < numberMasterColumns; i++) {
     6981                    double value = -elementAdd[start+i];
     6982                    if (fabs(value) > 1.0e-12) {
     6983                      infeas -= primal[i] * value;
     6984                      smallest = CoinMin(smallest, fabs(value));
     6985                      largest = CoinMax(largest, fabs(value));
     6986                      columnAdd[number] = i;
     6987                      elementAdd[number++] = value;
     6988                    }
     6989                  }
     6990                  // if elements large or small then scale?
     6991                  if (largest>1.0e8||smallest<1.0e-8) {
     6992                    sprintf(generalPrint,"For subproblem ray %d smallest - %g, largest %g - infeas %g",
     6993                           iBlock, smallest, largest, infeas);
     6994                    handler_->message(CLP_GENERAL2, messages_)
     6995                      << generalPrint
     6996                      << CoinMessageEol;
     6997                  }
     6998                  if (smallest<1.0e-12*largest) {
     6999                    sprintf(generalPrint,"Removing small elements");
     7000                    handler_->message(CLP_GENERAL2, messages_)
     7001                      << generalPrint
     7002                      << CoinMessageEol;
     7003                    smallest=1.0e-12*largest;
     7004                    int number2=number-1;
     7005                    number = start;
     7006                    for (int i = start; i < number2; i++) {
     7007                      double value = elementAdd[i];
     7008                      if (fabs(value) > smallest) {
     7009                        columnAdd[number] = columnAdd[i];
     7010                        elementAdd[number++] = value;
     7011                      }
     7012                    }
     7013                    columnAdd[number] = numberMasterColumns + iBlock;
     7014                    elementAdd[number++] = -1.0;
     7015                  }
     7016                  if (infeas < -1.0e-6||false) {
     7017                    double sum=0.0;
     7018                    for (int i=start;i<number;i++) {
     7019                      int iColumn=columnAdd[i];
     7020                      sum  += primal[iColumn]*elementAdd[i];
     7021                    }
     7022                    if (logLevel>2)
     7023                      printf("Sum %g rhs %g\n",sum,objValue);
     7024                    if (logLevel>1)
     7025                      printf("Cut for block %d has %d elements (infeasibility)\n",iBlock,number-start);
     7026                    blockPrint[numberProposals]=iBlock;
     7027                    // take
     7028                    objective[numberProposals] = objValue;
     7029                    rowAdd[++numberProposals] = number;
     7030                    when[numberRowsGenerated] = iPass;
     7031                    whichBlock[numberRowsGenerated++] = iBlock;
     7032                  }
     7033                } else {
     7034                  abort();
     7035                }
     7036              }
     7037            } else {
     7038              //printf("Can skip\n");
     7039            }
     7040            problemState[iBlock] &= ~4;
     7041          }
     7042          if (deletePrimal)
     7043            delete [] primal;
    53847044          if (numberProposals) {
    5385                master.addRows(numberProposals, NULL, objective,
     7045               sprintf(generalPrint,"%d cuts added with %d elements",
     7046                      numberProposals,rowAdd[numberProposals]);
     7047               handler_->message(CLP_GENERAL, messages_)
     7048                 << generalPrint
     7049                 << CoinMessageEol;
     7050               if (logLevel>2) {
     7051                 for (int i=0;i<numberProposals;i++) {
     7052                   printf("Cut %d block %d thetac %d ",i,blockPrint[i],
     7053                          blockPrint[i]+numberMasterColumns);
     7054                   int k=0;
     7055                   for (int j=rowAdd[i];j<rowAdd[i+1];j++) {
     7056                     if (k==12) {
     7057                       printf("\n");
     7058                       k=0;
     7059                     }
     7060                     k++;
     7061                     printf("(%d,%g) ",columnAdd[j],elementAdd[j]);
     7062                   }
     7063                   printf(" <= %g\n",objective[i]);
     7064                   //if (k)
     7065                   //printf("\n");
     7066                 }
     7067               }
     7068#ifdef TEST_MODEL
     7069               if (logLevel>3) {
     7070                 const double * solution = goodModel.primalColumnSolution();
     7071                 const double * solution2 = masterModel.primalColumnSolution();
     7072                 for (int i=0;i<numberProposals;i++) {
     7073                   double sum = 0.0;
     7074                   double sum2 = 0.0;
     7075                   for (int j=rowAdd[i];j<rowAdd[i+1];j++) {
     7076                     int iColumn=columnAdd[j];
     7077                     double value=elementAdd[j];
     7078                     sum += value*solution[iColumn];
     7079                     sum2 += value*solution2[iColumn];
     7080                   }
     7081                   if (sum2<objective[i]-1.0e-4) {
     7082                     sprintf(generalPrint,"Rhs for cut %d (from block %d) does not cutoff sum2 %g sum %g rhs %g)",
     7083                            i,blockPrint[i],sum2,sum,objective[i]);
     7084                     handler_->message(CLP_GENERAL2, messages_)
     7085                       << generalPrint
     7086                       << CoinMessageEol;
     7087                   }
     7088#define FIXUP_RHS 0
     7089                   if (sum>objective[i]+1.0e-4) {
     7090                     sprintf(generalPrint,"Rhs for cut %d (from block %d) is %g too low (rhs is %g)",
     7091                            i,blockPrint[i],sum-objective[i],objective[i]);
     7092                     handler_->message(CLP_GENERAL2, messages_)
     7093                       << generalPrint
     7094                       << CoinMessageEol;
     7095#if FIXUP_RHS == 1 || FIXUP_RHS ==3
     7096                     objective[i]=sum;
     7097#endif
     7098                   } else if (sum<objective[i]-1.0e-4) {
     7099                     sprintf(generalPrint,"Rhs for cut %d (from block %d) is %g ineffective (rhs is %g)",
     7100                            i,blockPrint[i],objective[i]-sum,objective[i]);
     7101                     handler_->message(CLP_GENERAL2, messages_)
     7102                       << generalPrint
     7103                       << CoinMessageEol;
     7104#if FIXUP_RHS == 2 || FIXUP_RHS ==3
     7105                     objective[i]=sum;
     7106#endif
     7107                   }
     7108                 }
     7109               }
     7110               if (logLevel>3) {
     7111                 goodModel.addRows(numberProposals, NULL, objective,
     7112                                   rowAdd, columnAdd, elementAdd);
     7113                 goodModel.dual();
     7114                 if (goodModel.problemStatus()==1||goodModel.objectiveValue()>goodValue+1.0e-5*fabs(goodValue)) {
     7115                   int numberRows=goodModel.numberRows();
     7116                   int numberStart=numberRows-numberProposals;
     7117                   double * upper = goodModel.rowUpper();
     7118                   for (int iRow=numberStart;iRow<numberRows;iRow++) {
     7119                     upper[iRow]=COIN_DBL_MAX;
     7120                   }
     7121                   for (int iRow=numberStart;iRow<numberRows;iRow++) {
     7122                     upper[iRow]=objective[iRow-numberStart];
     7123                     goodModel.allSlackBasis(true);
     7124                     goodModel.dual();
     7125                     if (goodModel.problemStatus()==1) {
     7126                       sprintf(generalPrint,"Cut %d makes infeasible - upper=%g",
     7127                              iRow-numberStart,upper[iRow]);
     7128                       handler_->message(CLP_GENERAL, messages_)
     7129                         << generalPrint
     7130                         << CoinMessageEol;
     7131                     } else if (goodModel.objectiveValue()>goodValue+1.0e-5*fabs(goodValue)) {
     7132                       sprintf(generalPrint,"Cut %d makes too expensive - upper=%g",
     7133                              iRow-numberStart,upper[iRow]);
     7134                       handler_->message(CLP_GENERAL, messages_)
     7135                         << generalPrint
     7136                         << CoinMessageEol;
     7137                       int iBlock=blockPrint[iRow-numberStart];
     7138                       ClpSimplex * temp = deBound(sub+iBlock);
     7139                       temp->allSlackBasis();
     7140                       temp->primal();
     7141                       // use ray
     7142                       double * solution = temp->infeasibilityRay();
     7143                       int numberRows2=temp->numberRows();
     7144                       // bug somewhere - if primal then flip
     7145                       if (true) {
     7146                         for (int i = 0; i < numberRows2; i++)
     7147                           solution[i] = - 1.0e-5 * solution[i];
     7148                       }
     7149                       double objValue=0.0;
     7150                       const double * lower = temp->rowLower();
     7151                       const double * upper = temp->rowUpper();
     7152                       for (int i = 0; i < numberRows2; i++) {
     7153                         if (solution[i] < -dualTolerance_) {
     7154                           // at upper
     7155                           assert (upper[i] < 1.0e30);
     7156                           objValue += solution[i] * upper[i];
     7157                         } else if (solution[i] > dualTolerance_) {
     7158                           // at lower
     7159                           assert (lower[i] > -1.0e30);
     7160                           objValue += solution[i] * lower[i];
     7161                         }
     7162                       }
     7163                       //printf("new objValue %g\n",objValue);
     7164                       int numberColumns2=sub[iBlock].numberColumns();
     7165                       double * temp2 = new double[numberColumns2];
     7166                       memset(temp2,0,numberColumns2*sizeof(double));
     7167                       sub[iBlock].clpMatrix()->transposeTimes(1.0,sub[iBlock].infeasibilityRay(), temp2);
     7168                       double loX=0.0;
     7169                       double upX=0.0;
     7170                       const double * lower2 = sub[iBlock].columnLower();
     7171                       const double * upper2 = sub[iBlock].columnUpper();
     7172                       for (int i=0;i<numberColumns2;i++) {
     7173                         if (sub[iBlock].getColumnStatus(i)!=basic&&
     7174                             fabs(sub[iBlock].primalColumnSolution()[i])>1.0e-5)
     7175                           printf("zz_inf %d has value %g\n",i,sub[iBlock].primalColumnSolution()[i]);
     7176                         if (sub[iBlock].getStatus(i)==atLowerBound) {
     7177                           loX += lower2[i]*temp2[i];
     7178                         } else if (sub[iBlock].getStatus(i)==atUpperBound) {
     7179                           upX += upper2[i]*temp2[i];
     7180                         }
     7181                       }
     7182                       printf("Offsets %g %g\n",loX,upX);
     7183                       memset(temp2,0,numberColumns2*sizeof(double));
     7184                       temp->clpMatrix()->transposeTimes(1.0,solution, temp2);
     7185                       loX=0.0;
     7186                       upX=0.0;
     7187                       lower2 = temp->columnLower();
     7188                       upper2 = temp->columnUpper();
     7189                       for (int i=0;i<numberColumns2;i++) {
     7190                         if (temp->getColumnStatus(i)!=basic&&
     7191                             fabs(temp->primalColumnSolution()[i])>1.0e-5)
     7192                           printf("zz_inf %d has value %g\n",i,temp->primalColumnSolution()[i]);
     7193                         if (temp->getStatus(i)==atLowerBound) {
     7194                           loX += lower2[i]*temp2[i];
     7195                         } else if (temp->getStatus(i)==atUpperBound) {
     7196                           upX += upper2[i]*temp2[i];
     7197                         }
     7198                       }
     7199                       printf("Offsets %g %g\n",loX,upX);
     7200                       delete [] temp2;
     7201                       delete temp;
     7202                     }
     7203                     upper[iRow]=COIN_DBL_MAX;
     7204                   }
     7205                 }
     7206                 double objValue=goodModel.objectiveValue();
     7207                 const double * obj = goodModel.objective();
     7208                 const double * solution = goodModel.primalColumnSolution();
     7209                 double obj1=0.0;
     7210                 for (int i=0;i<numberMasterColumns;i++)
     7211                   obj1 += obj[i]*solution[i];
     7212                 double obj2=0.0;
     7213                 for (int i=numberMasterColumns;i<numberMasterColumns+numberBlocks;i++)
     7214                   obj2 += obj[i]*solution[i];
     7215                 double obj3=0.0;
     7216                 for (int i=numberMasterColumns+numberBlocks;i<goodModel.numberColumns();i++)
     7217                   obj3 += obj[i]*solution[i];
     7218                 //assert (fabs(goodValue-objValue)<1.0e-3+1.0e-7*fabs(goodValue));
     7219                 printf("XXXX good %g this %g difference %g - objs %g, %g, %g\n",
     7220                        goodValue,objValue,goodValue-objValue,obj1,obj2,obj3);
     7221               }
     7222#endif
     7223               masterModel.addRows(numberProposals, NULL, objective,
    53867224                              rowAdd, columnAdd, elementAdd);
    5387           }
    5388      }
    5389      std::cout << "Time at end of Benders " << CoinCpuTime() - time1 << " seconds" << std::endl;
    5390      //master.scaling(0);
    5391      //master.primal(1);
    5392      loadProblem(*model);
     7225          }
     7226     }
     7227     sprintf(generalPrint,"Time at end of Benders %.2f seconds",CoinCpuTime() - time1);
     7228     handler_->message(CLP_GENERAL, messages_)
     7229       << generalPrint
     7230       << CoinMessageEol;
     7231     delete [] problemState;
     7232     for (int iBlock = 0; iBlock < numberBlocks; iBlock++) {
     7233       delete [] modification[iBlock];
     7234     }
     7235     delete [] modification;
     7236#ifdef ADD_ARTIFICIALS
     7237     delete [] originalSubColumns;
     7238     for (int i=0;i<numberBlocks;i++)
     7239       delete [] saveObjective[i];
     7240     delete [] saveObjective;
     7241#endif
     7242     //masterModel.scaling(0);
     7243     //masterModel.primal(1);
     7244     if (!options.independentOption(1))
     7245       loadProblem(*model);
    53937246     // now put back a good solution
    5394      const double * columnSolution = master.primalColumnSolution();
     7247     const double * columnSolution = masterModel.primalColumnSolution();
    53957248     double * fullColumnSolution = primalColumnSolution();
    53967249     const int * columnBack = model->coinBlock(masterBlock)->originalColumns();
     
    53987251     const int * rowBack = model->coinBlock(masterBlock)->originalRows();
    53997252     int numberRows2 = model->coinBlock(masterBlock)->numberRows();
     7253#ifndef NDEBUG
     7254     for (int iColumn = 0; iColumn < numberColumns_; iColumn++)
     7255       fullColumnSolution[iColumn]=COIN_DBL_MAX;
     7256#endif
    54007257     for (int iColumn = 0; iColumn < numberColumns2; iColumn++) {
    54017258          int kColumn = columnBack[iColumn];
    5402           setColumnStatus(kColumn, master.getColumnStatus(iColumn));
     7259          setColumnStatus(kColumn, masterModel.getColumnStatus(iColumn));
    54037260          fullColumnSolution[kColumn] = columnSolution[iColumn];
    54047261     }
    54057262     for (int iRow = 0; iRow < numberRows2; iRow++) {
    54067263          int kRow = rowBack[iRow];
    5407           setStatus(kRow, master.getStatus(iRow));
     7264          setRowStatus(kRow, masterModel.getRowStatus(iRow));
    54087265          //fullSolution[kRow]=solution[iRow];
    54097266     }
     
    54197276               int kColumn = columnBack[iColumn];
    54207277               setColumnStatus(kColumn, sub[iBlock].getColumnStatus(iColumn));
     7278#ifndef NDEBUG
     7279               assert(fullColumnSolution[kColumn]==COIN_DBL_MAX);
     7280#endif
    54217281               fullColumnSolution[kColumn] = subColumnSolution[iColumn];
    54227282          }
    54237283          for (int iRow = 0; iRow < numberRows2; iRow++) {
    54247284               int kRow = rowBack[iRow];
    5425                setStatus(kRow, sub[iBlock].getStatus(iRow));
    5426                setStatus(kRow, atLowerBound);
    5427           }
    5428      }
     7285               setRowStatus(kRow, sub[iBlock].getRowStatus(iRow));
     7286               //setStatus(kRow, atLowerBound);
     7287          }
     7288     }
     7289#ifndef NDEBUG
     7290     for (int iColumn = 0; iColumn < numberColumns_; iColumn++)
     7291       assert(fullColumnSolution[iColumn]!=COIN_DBL_MAX);
     7292#endif
    54297293     double * fullSolution = primalRowSolution();
    54307294     CoinZeroN(fullSolution, numberRows_);
    54317295     times(1.0, fullColumnSolution, fullSolution);
    5432      //int numberColumns=model->numberColumns();
    5433      //for (int iColumn=0;iColumn<numberColumns;iColumn++)
    5434      //setColumnStatus(iColumn,ClpSimplex::superBasic);
    5435      std::cout << "Time before cleanup of full model " << CoinCpuTime() - time1 << " seconds" << std::endl;
     7296     int numberRowBasic=0;
     7297#ifndef NDEBUG
     7298     int numberInfeasibilities=0;
     7299     double sumInfeasibilities=0.0;
     7300#endif
     7301     for (int iRow=0;iRow<numberRows_;iRow++) {
     7302       if (getRowStatus(iRow)==ClpSimplex::basic)
     7303         numberRowBasic++;
     7304#ifndef NDEBUG
     7305       if (fullSolution[iRow]<rowLower_[iRow]-primalTolerance_) {
     7306         numberInfeasibilities++;
     7307         sumInfeasibilities -= fullSolution[iRow]-rowLower_[iRow];
     7308         if (getRowStatus(iRow)!=basic)
     7309           setRowStatus(iRow,superBasic);
     7310       } else if (fullSolution[iRow]>rowUpper_[iRow]+primalTolerance_) {
     7311         numberInfeasibilities++;
     7312         sumInfeasibilities += fullSolution[iRow]-rowUpper_[iRow];
     7313         if (getRowStatus(iRow)!=basic)
     7314           setRowStatus(iRow,superBasic);
     7315       }
     7316#endif
     7317     }
     7318     int numberColumnBasic=0;
     7319     for (int iColumn=0;iColumn<numberColumns_;iColumn++)
     7320       if (getColumnStatus(iColumn)==ClpSimplex::basic)
     7321         numberColumnBasic++;
     7322     sprintf(generalPrint,"%d row basic %d col basic (total %d) - wanted %d",
     7323            numberRowBasic,numberColumnBasic,
     7324            numberRowBasic+numberColumnBasic,numberRows_);
     7325     handler_->message(CLP_GENERAL2, messages_)
     7326       << generalPrint
     7327       << CoinMessageEol;
     7328#ifndef NDEBUG
     7329     sprintf(generalPrint,"%d infeasibilities summing to %g",
     7330            numberInfeasibilities,sumInfeasibilities);
     7331     handler_->message(CLP_GENERAL2, messages_)
     7332       << generalPrint
     7333       << CoinMessageEol;
     7334#endif
     7335     //for (int i=0;i<numberRows_;i++) setRowStatus(i,basic);
     7336     sprintf(generalPrint,"Time before cleanup of full model %.2f seconds",CoinCpuTime() - time1);
     7337     handler_->message(CLP_GENERAL, messages_)
     7338       << generalPrint
     7339       << CoinMessageEol;
     7340#ifdef ABC_INHERIT
     7341     this->dealWithAbc(1,1,true);
     7342#else
    54367343     this->primal(1);
    5437      std::cout << "Total time " << CoinCpuTime() - time1 << " seconds" << std::endl;
     7344#endif
     7345     sprintf(generalPrint,"Total time %.2f seconds",CoinCpuTime() - time1);
     7346     handler_->message(CLP_GENERAL, messages_)
     7347       << generalPrint
     7348       << CoinMessageEol;
    54387349     delete [] rowCounts;
    54397350     //delete [] sol;
  • trunk/Clp/src/ClpSolve.hpp

    r1924 r2025  
    3030          useBarrierNoCross,
    3131          automatic,
     32          tryDantzigWolfe,
     33          tryBenders,
    3234          notImplemented
    3335     };
     
    222224          independentOptions_[2] = value;
    223225     }
     226     inline void setIndependentOption(int type,int value) {
     227          independentOptions_[type]  = value;
     228     }
     229     inline int independentOption(int type) const {
     230          return independentOptions_[type];
     231     }
    224232     //@}
    225233
     
    244252         1 - To be copied over to presolve options
    245253         2 - max substitution level
     254         If Dantzig Wolfe/benders 0 is number blocks, 2 is #passes (notional)
    246255     */
    247256     int independentOptions_[3];
     
    294303     /// Returns real primal infeasibility (if -1) - current if (0)
    295304     double lastInfeasibility(int back = 1) const;
     305     /// Returns number of primal infeasibilities (if -1) - current if (0)
     306     int numberInfeasibilities(int back = 1) const;
    296307     /// Modify objective e.g. if dual infeasible in dual
    297308     void modifyObjective(double value);
  • trunk/Clp/src/OsiClp/OsiClpSolverInterface.cpp

    r1972 r2025  
    22382238      }
    22392239    }
     2240    //#define REDO_STEEPEST_EDGE
     2241#ifdef REDO_STEEPEST_EDGE
     2242    if (dynamic_cast<ClpDualRowSteepest *>(modelPtr_->dualRowPivot())) {
     2243      // new copy (should think about keeping around)
     2244      ClpDualRowSteepest steepest;
     2245      modelPtr_->setDualRowPivotAlgorithm(steepest);
     2246    }     
     2247#endif
    22402248    // Start of fast iterations
    22412249    bool alwaysFinish= ((specialOptions_&32)==0) ? true : false;
     
    24432451    int status;
    24442452    if (!infeasible) {
     2453#ifdef REDO_STEEPEST_EDGE
     2454      if (dynamic_cast<ClpDualRowSteepest *>(smallModel_->dualRowPivot())) {
     2455        // new copy
     2456        ClpDualRowSteepest steepest;
     2457        smallModel_->setDualRowPivotAlgorithm(steepest);
     2458      }     
     2459#endif
    24452460      status = static_cast<ClpSimplexDual *>(smallModel_)->fastDual(alwaysFinish);
    24462461    } else {
Note: See TracChangeset for help on using the changeset viewer.