Changeset 854 for branches


Ignore:
Timestamp:
Sep 15, 2006 4:48:56 PM (13 years ago)
Author:
forrest
Message:

many changes

Location:
branches/devel/Clp/src
Files:
26 edited

Legend:

Unmodified
Added
Removed
  • branches/devel/Clp/src/CbcOrClpParam.cpp

    r837 r854  
    764764      model.setIntParam(CbcModel::CbcMaxNumNode,value);
    765765      break;
     766    case MAXSOLS:
     767      oldValue=model.getIntParam(CbcModel::CbcMaxNumSol);
     768      model.setIntParam(CbcModel::CbcMaxNumSol,value);
     769      break;
    766770    case STRONGBRANCHING:
    767771      oldValue=model.numberStrong();
     
    802806  case MAXNODES:
    803807    value = model.getIntParam(CbcModel::CbcMaxNumNode);
     808    break;
     809  case MAXSOLS:
     810    value = model.getIntParam(CbcModel::CbcMaxNumSol);
    804811    break;
    805812  case STRONGBRANCHING:
     
    11601167First just try with default settings and look carefully at the log file.  Did cuts help?  Did they take too long?  \
    11611168Look at output to see which cuts were effective and then do some tuning.  You will see that the \
    1162 options for cuts are off, on, root and ifmove.  Off is \
     1169options for cuts are off, on, root and ifmove, forceon.  Off is \
    11631170obvious, on means that this cut generator will be tried in the branch and cut tree (you can fine tune using \
    11641171'depth').  Root means just at the root node while 'ifmove' means that cuts will be used in the tree if they \
    1165 look as if they are doing some good and moving the objective value.  If pre-processing reduced the size of the \
     1172look as if they are doing some good and moving the objective value.  Forceon is same as on but forces code to use \
     1173cut generator at every node.  For probing forceonbut just does fixing probing in tree - not strengthening etc.  \
     1174If pre-processing reduced the size of the \
    11661175problem or strengthened many coefficients then it is probably wise to leave it on.  Switch off heuristics \
    11671176which did not provide solutions.  The other major area to look at is the search.  Hopefully good solutions \
     
    12141223  parameters[numberParameters-1].append("root");
    12151224  parameters[numberParameters-1].append("ifmove");
     1225  parameters[numberParameters-1].append("forceOn");
    12161226  parameters[numberParameters-1].setLonghelp
    12171227    (
     
    12271237     "This switches on a heuristic which does branch and cut on the problem given by just \
    12281238using variables which have appeared in one or more solutions. \
    1229 It is obviously only tries after two or more solutions."
     1239It obviously only tries after two or more solutions."
    12301240     );
    12311241  parameters[numberParameters++]=
     
    12341244  parameters[numberParameters-1].append("pri!orities");
    12351245  parameters[numberParameters-1].append("column!Order?");
     1246  parameters[numberParameters-1].append("01f!irst?");
     1247  parameters[numberParameters-1].append("01l!ast?");
     1248  parameters[numberParameters-1].append("length!?");
    12361249  parameters[numberParameters-1].setLonghelp
    12371250    (
     
    13111324  parameters[numberParameters-1].append("root");
    13121325  parameters[numberParameters-1].append("ifmove");
     1326  parameters[numberParameters-1].append("forceOn");
    13131327  parameters[numberParameters-1].setLonghelp
    13141328    (
     
    14851499    parameters[numberParameters-1].append("root");
    14861500    parameters[numberParameters-1].append("ifmove");
    1487   parameters[numberParameters-1].setLonghelp
     1501    parameters[numberParameters-1].append("forceOn");
     1502    parameters[numberParameters-1].setLonghelp
    14881503    (
    14891504     "This switches on flow cover cuts (either at root or in entire tree) \
     
    15161531  parameters[numberParameters-1].append("root");
    15171532  parameters[numberParameters-1].append("ifmove");
     1533  parameters[numberParameters-1].append("forceOn");
    15181534  parameters[numberParameters-1].setLonghelp
    15191535    (
     
    16391655  parameters[numberParameters-1].append("root");
    16401656  parameters[numberParameters-1].append("ifmove");
     1657  parameters[numberParameters-1].append("forceOn");
    16411658  parameters[numberParameters-1].setLonghelp
    16421659    (
    16431660     "This switches on knapsack cuts (either at root or in entire tree) \
     1661See branchAndCut for information on options."
     1662     );
     1663  parameters[numberParameters++]=
     1664    CbcOrClpParam("lift!AndProjectCuts","Whether to use Lift and Project cuts",
     1665                  "off",LANDPCUTS);
     1666  parameters[numberParameters-1].append("on");
     1667  parameters[numberParameters-1].append("root");
     1668  parameters[numberParameters-1].append("ifmove");
     1669  parameters[numberParameters-1].append("forceOn");
     1670  parameters[numberParameters-1].setLonghelp
     1671    (
     1672     "Lift and project cuts - may be expensive to compute. \
    16441673See branchAndCut for information on options."
    16451674     );
     
    17091738     "This is a repeatable way to limit search.  Normally using time is easier \
    17101739but then the results may not be repeatable."
     1740     );
     1741  parameters[numberParameters++]=
     1742    CbcOrClpParam("maxS!olutions","Maximum number of solutions to get",
     1743                  1,2147483647,MAXSOLS);
     1744  parameters[numberParameters-1].setLonghelp
     1745    (
     1746     "You may want to stop after (say) two solutions or an hour."
    17111747     );
    17121748#endif
     
    17331769  parameters[numberParameters-1].append("root");
    17341770  parameters[numberParameters-1].append("ifmove");
     1771  parameters[numberParameters-1].append("forceOn");
    17351772  parameters[numberParameters-1].setLonghelp
    17361773    (
     
    19491986  parameters[numberParameters-1].append("sos");
    19501987  parameters[numberParameters-1].append("trysos");
     1988  parameters[numberParameters-1].append("equalall");
    19511989  parameters[numberParameters-1].append("strategy");
    19521990  parameters[numberParameters-1].setLonghelp
     
    19561994 Save option saves on file presolved.mps.  equal will turn <= cliques into \
    19571995==.  sos will create sos sets if all 0-1 in sets (well one extra is allowed) \
    1958 and no overlaps.  trysos is same but allows any number extra.  strategy is as \
     1996and no overlaps.  trysos is same but allows any number extra.  equalall will turn all \
     1997valid inequalities into equalities with integer slacks.  strategy is as \
    19591998on but uses CbcStrategy."
    19601999     );
     
    20682107  parameters[numberParameters-1].append("root");
    20692108  parameters[numberParameters-1].append("ifmove");
     2109  parameters[numberParameters-1].append("forceOn");
     2110  parameters[numberParameters-1].append("forceOnBut");
    20702111  parameters[numberParameters-1].setLonghelp
    20712112    (
     
    21142155    parameters[numberParameters-1].append("root");
    21152156    parameters[numberParameters-1].append("ifmove");
    2116   parameters[numberParameters-1].setLonghelp
     2157    parameters[numberParameters-1].append("forceOn");
     2158    parameters[numberParameters-1].setLonghelp
    21172159    (
    21182160     "This switches on reduce and split  cuts (either at root or in entire tree) \
     
    23732415#ifdef COIN_HAS_CBC
    23742416  parameters[numberParameters++]=
     2417    CbcOrClpParam("tune!PreProcess","Dubious tuning parameters",
     2418                  0,1000000,PROCESSTUNE,false);
     2419  parameters[numberParameters-1].setLonghelp
     2420    (
     2421     "For making equality cliques this is minimumsize.  Also for adding \
     2422integer slacks.  May be used for more later"
     2423     );
     2424  parameters[numberParameters++]=
    23752425    CbcOrClpParam("two!MirCuts","Whether to use Two phase Mixed Integer Rounding cuts",
    23762426                  "off",TWOMIRCUTS);
     
    23782428  parameters[numberParameters-1].append("root");
    23792429  parameters[numberParameters-1].append("ifmove");
     2430  parameters[numberParameters-1].append("forceOn");
    23802431  parameters[numberParameters-1].setLonghelp
    23812432    (
     
    24092460     );
    24102461#endif
     2462  parameters[numberParameters++]=
     2463    CbcOrClpParam("vector","Whether to use vector? Form of matrix in simplex",
     2464                  "off",VECTOR,7,false);
     2465  parameters[numberParameters-1].append("on");
     2466  parameters[numberParameters-1].setLonghelp
     2467    (
     2468     "If this is on and ClpPackedMatrix uses extra column copy in odd format."
     2469     );
    24112470  parameters[numberParameters++]=
    24122471    CbcOrClpParam("verbose","Switches on longer help on single ?",
  • branches/devel/Clp/src/CbcOrClpParam.hpp

    r800 r854  
    6060    MAXFACTOR,PERTVALUE,MAXITERATION,PRESOLVEPASS,IDIOT,SPRINT,
    6161    OUTPUTFORMAT,SLPVALUE,PRESOLVEOPTIONS,PRINTOPTIONS,SPECIALOPTIONS,
    62     SUBSTITUTION,DUALIZE,VERBOSE,THREADS,CPP,
     62    SUBSTITUTION,DUALIZE,VERBOSE,THREADS,CPP,CUTPASS,PROCESSTUNE,
    6363
    6464    STRONGBRANCHING=151,CUTDEPTH, MAXNODES,NUMBERBEFORE,NUMBERANALYZE,
    65     NUMBERMINI,MIPOPTIONS,MOREMIPOPTIONS,MAXHOTITS,FPUMPITS,CUTPASS,
     65    NUMBERMINI,MIPOPTIONS,MOREMIPOPTIONS,MAXHOTITS,FPUMPITS,MAXSOLS,
    6666#ifdef COIN_HAS_CBC
    6767    LOGLEVEL ,
     
    7070    DIRECTION=201,DUALPIVOT,SCALING,ERRORSALLOWED,KEEPNAMES,SPARSEFACTOR,
    7171    PRIMALPIVOT,PRESOLVE,CRASH,BIASLU,PERTURBATION,MESSAGES,AUTOSCALE,
    72     CHOLESKY,KKT,BARRIERSCALE,GAMMA,CROSSOVER,PFI,INTPRINT,
     72    CHOLESKY,KKT,BARRIERSCALE,GAMMA,CROSSOVER,PFI,INTPRINT,VECTOR,
    7373   
    7474    NODESTRATEGY = 251,BRANCHSTRATEGY,CUTSSTRATEGY,HEURISTICSTRATEGY,
     
    7676    ROUNDING,SOLVER,CLIQUECUTS,COSTSTRATEGY,FLOWCUTS,MIXEDCUTS,
    7777    TWOMIRCUTS,PREPROCESS,FPUMP,GREEDY,COMBINE,LOCALTREE,USESOLUTION,SOS,
     78    LANDPCUTS,
    7879   
    7980    DIRECTORY=301,IMPORT,EXPORT,RESTORE,SAVE,DUALSIMPLEX,PRIMALSIMPLEX,EITHERSIMPLEX,
  • branches/devel/Clp/src/ClpDynamicMatrix.cpp

    r800 r854  
    230230  numberElements_ = CoinMin(numberElements_,numberElements)+originalMatrix->getNumElements();
    231231  matrix_ = originalMatrix;
    232   zeroElements_ = false;
     232  flags_ &= ~1;
    233233  // resize model (matrix stays same)
    234234  int newRowSize = numberRows+CoinMin(numberSets_,CoinMax(frequency,numberRows))+1;
  • branches/devel/Clp/src/ClpEventHandler.cpp

    r754 r854  
    6767  model_= model;
    6868}
     69//#############################################################################
     70// Constructors / Destructor / Assignment
     71//#############################################################################
     72
     73//-------------------------------------------------------------------
     74// Default Constructor
     75//-------------------------------------------------------------------
     76ClpDisasterHandler::ClpDisasterHandler (ClpSimplex * model) :
     77  model_(model)
     78{
     79
     80}
     81
     82//-------------------------------------------------------------------
     83// Copy constructor
     84//-------------------------------------------------------------------
     85ClpDisasterHandler::ClpDisasterHandler (const ClpDisasterHandler & rhs)
     86  : model_(rhs.model_)
     87
     88}
     89
     90//-------------------------------------------------------------------
     91// Destructor
     92//-------------------------------------------------------------------
     93ClpDisasterHandler::~ClpDisasterHandler ()
     94{
     95}
     96
     97//----------------------------------------------------------------
     98// Assignment operator
     99//-------------------------------------------------------------------
     100ClpDisasterHandler &
     101ClpDisasterHandler::operator=(const ClpDisasterHandler& rhs)
     102{
     103  if (this != &rhs) {
     104    model_ = rhs.model_;
     105  }
     106  return *this;
     107}
     108/* set model. */
     109void
     110ClpDisasterHandler::setSimplex(ClpSimplex * model)
     111{
     112  model_= model;
     113}
    69114
    70115
  • branches/devel/Clp/src/ClpEventHandler.hpp

    r754 r854  
    8686  //@}
    8787};
     88/** Base class for Clp disaster handling
     89   
     90This is here to allow for disaster handling.  By disaster I mean that Clp
     91would otherwise give up
     92
     93*/
     94
     95class ClpDisasterHandler  {
     96 
     97public:
     98  /**@name Virtual methods that the derived classe should provide.
     99  */
     100  //@{
     101  /// Into simplex
     102  virtual void intoSimplex()=0;
     103  /// Checks if disaster
     104  virtual bool check() const = 0;
     105  /// saves information for next attempt
     106  virtual void saveInfo() =0;
     107  //@}
     108 
     109 
     110  /**@name Constructors, destructor */
     111
     112  //@{
     113  /** Default constructor. */
     114  ClpDisasterHandler(ClpSimplex * model = NULL);
     115  /** Destructor */
     116  virtual ~ClpDisasterHandler();
     117  // Copy
     118  ClpDisasterHandler(const ClpDisasterHandler&);
     119  // Assignment
     120  ClpDisasterHandler& operator=(const ClpDisasterHandler&);
     121  /// Clone
     122  virtual ClpDisasterHandler * clone() const =0;
     123
     124  //@}
     125 
     126  /**@name Sets/gets */
     127
     128  //@{
     129  /** set model. */
     130  void setSimplex(ClpSimplex * model);
     131  /// Get model
     132  inline ClpSimplex * simplex() const
     133  { return model_;};
     134  //@}
     135 
     136 
     137protected:
     138  /**@name Data members
     139     The data members are protected to allow access for derived classes. */
     140  //@{
     141  /// Pointer to simplex
     142  ClpSimplex * model_;
     143  //@}
     144};
    88145#endif
  • branches/devel/Clp/src/ClpFactorization.cpp

    r754 r854  
    247247      }
    248248      if (numberBasic>model->maximumBasic()) {
    249 #ifndef NDEBUG
     249#if 0 // ndef NDEBUG
    250250        printf("%d basic - should only be %d\n",
    251251               numberBasic,numberRows);
  • branches/devel/Clp/src/ClpGubDynamicMatrix.cpp

    r800 r854  
    168168  numberElements_ = CoinMin(numberElements_,numberElements)+originalMatrix->getNumElements();
    169169  matrix_ = originalMatrix;
    170   zeroElements_ = false;
     170  flags_ &= ~1;
    171171  // resize model (matrix stays same)
    172172  model->resize(numberRows,numberNeeded);
  • branches/devel/Clp/src/ClpGubMatrix.cpp

    r800 r854  
    40734073// Correct sequence in and out to give true value
    40744074void
    4075 ClpGubMatrix::correctSequence(int & sequenceIn, int & sequenceOut) const
     4075ClpGubMatrix::correctSequence(const ClpSimplex * model,int & sequenceIn, int & sequenceOut)
    40764076{
    4077   sequenceIn = trueSequenceIn_;
    4078   sequenceOut = trueSequenceOut_;
     4077  if (sequenceIn!=-999) {
     4078    sequenceIn = trueSequenceIn_;
     4079    sequenceOut = trueSequenceOut_;
     4080  }
    40794081}
  • branches/devel/Clp/src/ClpGubMatrix.hpp

    r754 r854  
    153153  virtual int synchronize(ClpSimplex * model,int mode);
    154154  /// Correct sequence in and out to give true value
    155   virtual void correctSequence(int & sequenceIn, int & sequenceOut) const;
     155  virtual void correctSequence(const ClpSimplex * model,int & sequenceIn, int & sequenceOut) ;
    156156  //@}
    157157
  • branches/devel/Clp/src/ClpMain.cpp

    r828 r854  
    6060
    6161int mainTest (int argc, const char *argv[],int algorithm,
    62               ClpSimplex empty, bool doPresolve,int switchOff);
     62              ClpSimplex empty, bool doPresolve,int switchOff,bool doVector);
    6363static void statistics(ClpSimplex * originalModel, ClpSimplex * model);
    6464static void generateCode(const char * fileName,int type);
     
    8686    int presolveOptions=0;
    8787    int doCrash=0;
     88    int doVector=0;
    8889    int doSprint=-1;
    8990    // set reasonable defaults
     
    475476              doCrash=action;
    476477              break;
     478            case VECTOR:
     479              doVector=action;
     480              break;
    477481            case MESSAGES:
    478482              models[iModel].messageHandler()->setPrefix(action!=0);
     
    590594              solveOptions.setSpecialOption(4,presolveOptions);
    591595              solveOptions.setSpecialOption(5,printOptions&1);
     596              if (doVector) {
     597                ClpMatrixBase * matrix = models[iModel].clpMatrix();
     598                if (dynamic_cast< ClpPackedMatrix*>(matrix)) {
     599                  ClpPackedMatrix * clpMatrix = dynamic_cast< ClpPackedMatrix*>(matrix);
     600                  clpMatrix->makeSpecialColumnCopy();
     601                }
     602              }
    592603              if (method==ClpSolve::useDual) {
    593604                // dual
     
    13791390              models[iModel].setSpecialOptions(0);
    13801391              mainTest(nFields,fields,algorithm,models[iModel],
    1381                        (preSolve!=0),specialOptions);
     1392                       (preSolve!=0),specialOptions,doVector!=0);
    13821393            }
    13831394            break;
     
    13981409              if (models[iModel].numberRows())
    13991410                algorithm=7;
    1400               mainTest(nFields,fields,algorithm,models[iModel],(preSolve!=0),specialOptions);
     1411              mainTest(nFields,fields,algorithm,models[iModel],(preSolve!=0),specialOptions,doVector!=0);
    14011412            }
    14021413            break;
  • branches/devel/Clp/src/ClpMatrixBase.cpp

    r800 r854  
    548548// Correct sequence in and out to give true value
    549549void
    550 ClpMatrixBase::correctSequence(int & sequenceIn, int & sequenceOut) const
     550ClpMatrixBase::correctSequence(const ClpSimplex * model,int & sequenceIn, int & sequenceOut)
    551551{
    552552}
  • branches/devel/Clp/src/ClpMatrixBase.hpp

    r800 r854  
    227227      mode=12  - after factorize but before permute stuff
    228228      mode=13  - at end of simplex to delete stuff
     229
    229230  */
    230231  virtual int generalExpanded(ClpSimplex * model,int mode,int & number);
     
    242243  /// Returns reduced cost of a variable
    243244  double reducedCost(ClpSimplex * model,int sequence) const;
    244   /// Correct sequence in and out to give true value
    245   virtual void correctSequence(int & sequenceIn, int & sequenceOut) const;
     245  /// Correct sequence in and out to give true value (if both -1 maybe do whole matrix)
     246  virtual void correctSequence(const ClpSimplex * model,int & sequenceIn, int & sequenceOut) ;
    246247  //@}
    247248 
  • branches/devel/Clp/src/ClpModel.cpp

    r847 r854  
    30013001  int iRow;
    30023002  for (iRow=first; iRow<last;iRow++) {
    3003     rowNames_[iRow]= rowNames[iRow-first];
    3004     maxLength = CoinMax(maxLength,(unsigned int) strlen(rowNames[iRow-first]));
     3003    if (rowNames[iRow-first]&&strlen(rowNames[iRow-first])) {
     3004      rowNames_[iRow]= rowNames[iRow-first];
     3005      maxLength = CoinMax(maxLength,(unsigned int) strlen(rowNames[iRow-first]));
     3006    } else {
     3007      maxLength = CoinMax(maxLength,(unsigned int) 8);
     3008      char name[9];
     3009      sprintf(name,"R%7.7d",iRow);
     3010      rowNames_[iRow]=name;
     3011    }
    30053012  }
    30063013  // May be too big - but we would have to check both rows and columns to be exact
     
    30173024  int iColumn;
    30183025  for (iColumn=first; iColumn<last;iColumn++) {
    3019     columnNames_[iColumn]= columnNames[iColumn-first];
    3020     maxLength = CoinMax(maxLength,(unsigned int) strlen(columnNames[iColumn-first]));
     3026    if (columnNames[iColumn-first]&&strlen(columnNames[iColumn-first])) {
     3027      columnNames_[iColumn]= columnNames[iColumn-first];
     3028      maxLength = CoinMax(maxLength,(unsigned int) strlen(columnNames[iColumn-first]));
     3029    } else {
     3030      maxLength = CoinMax(maxLength,(unsigned int) 8);
     3031      char name[9];
     3032      sprintf(name,"C%7.7d",iColumn);
     3033      columnNames_[iColumn]=name;
     3034    }
    30213035  }
    30223036  // May be too big - but we would have to check both rows and columns to be exact
  • branches/devel/Clp/src/ClpPackedMatrix.cpp

    r839 r854  
    88#include "CoinIndexedVector.hpp"
    99#include "CoinHelperFunctions.hpp"
     10//#define THREAD
    1011
    1112#include "ClpSimplex.hpp"
     
    1516#include "ClpQuadraticObjective.hpp"
    1617#endif
    17 //#define THREAD
    1818// at end to get min/max!
    1919#include "ClpPackedMatrix.hpp"
     
    4444    matrix_(NULL),
    4545    numberActiveColumns_(0),
    46     zeroElements_(false),
    47     hasGaps_(true),
    48     rowCopy_(NULL)
     46    flags_(2),
     47    rowCopy_(NULL),
     48    columnCopy_(NULL)
    4949{
    5050  setType(1);
     
    5959  matrix_ = new CoinPackedMatrix(*(rhs.matrix_));
    6060  numberActiveColumns_ = rhs.numberActiveColumns_;
    61   zeroElements_ = rhs.zeroElements_;
    62   hasGaps_ = rhs.hasGaps_;
     61  flags_ = rhs.flags_;
    6362  int numberRows = getNumRows();
    6463  if (rhs.rhsOffset_&&numberRows) {
     
    6867  }
    6968  if (rhs.rowCopy_) {
     69    assert ((flags_&4)!=0);
    7070    rowCopy_ = new ClpPackedMatrix2(*rhs.rowCopy_);
    7171  } else {
    7272    rowCopy_ = NULL;
     73  }
     74  if (rhs.columnCopy_) {
     75    assert ((flags_&8)!=0);
     76    columnCopy_ = new ClpPackedMatrix3(*rhs.columnCopy_);
     77  } else {
     78    columnCopy_ = NULL;
    7379  }
    7480}
     
    8187
    8288  matrix_ = rhs;
    83   zeroElements_ = false;
    84   hasGaps_ = true;
     89  flags_ = 2;
    8590  numberActiveColumns_ = matrix_->getNumCols();
    8691  rowCopy_ = NULL;
     92  columnCopy_=NULL;
    8793  setType(1);
    8894 
     
    95101  numberActiveColumns_ = matrix_->getNumCols();
    96102  rowCopy_ = NULL;
    97   zeroElements_ = false;
    98   hasGaps_ = true;
     103  flags_ = 2;
     104  columnCopy_=NULL;
    99105  setType(1);
    100106 
     
    108114  delete matrix_;
    109115  delete rowCopy_;
     116  delete columnCopy_;
    110117}
    111118
     
    121128    matrix_ = new CoinPackedMatrix(*(rhs.matrix_));
    122129    numberActiveColumns_ = rhs.numberActiveColumns_;
    123     zeroElements_ = rhs.zeroElements_;
    124     hasGaps_ = rhs.hasGaps_;
     130    flags_ = rhs.flags_;
    125131    delete rowCopy_;
     132    delete columnCopy_;
    126133    if (rhs.rowCopy_) {
     134      assert ((flags_&4)!=0);
    127135      rowCopy_ = new ClpPackedMatrix2(*rhs.rowCopy_);
    128136    } else {
    129137      rowCopy_ = NULL;
     138    }
     139    if (rhs.columnCopy_) {
     140      assert ((flags_&8)!=0);
     141      columnCopy_ = new ClpPackedMatrix3(*rhs.columnCopy_);
     142    } else {
     143      columnCopy_ = NULL;
    130144    }
    131145  }
     
    160174                                 numberColumns,whichColumns);
    161175  numberActiveColumns_ = matrix_->getNumCols();
    162   zeroElements_ = rhs.zeroElements_;
    163   hasGaps_ = rhs.hasGaps_;
    164176  rowCopy_ = NULL;
     177  flags_ = rhs.flags_;
     178  columnCopy_=NULL;
    165179}
    166180ClpPackedMatrix::ClpPackedMatrix (
     
    173187                                 numberColumns,whichColumns);
    174188  numberActiveColumns_ = matrix_->getNumCols();
    175   zeroElements_ = false;
    176   hasGaps_=true;
    177189  rowCopy_ = NULL;
     190  flags_ = 2;
     191  columnCopy_=NULL;
    178192  setType(1);
    179193}
     
    190204  //copy->matrix_->removeGaps();
    191205  copy->numberActiveColumns_ = copy->matrix_->getNumCols();
    192   copy->hasGaps_=false;
     206  copy->flags_ = flags_&(~2); // no gaps
    193207  return copy;
    194208}
     
    227241  const int * columnLength = matrix_->getVectorLengths();
    228242  const double * elementByColumn = matrix_->getElements();
    229   if (!hasGaps_) {
     243  if (!(flags_&2)) {
    230244    if (scalar==1.0) {
    231245      CoinBigIndex start=columnStart[0];
     
    326340    const double * elementByColumn = matrix_->getElements();
    327341    if (!spare) {
    328       if (!hasGaps_) {
     342      if (!(flags_&2)) {
    329343        CoinBigIndex start=columnStart[0];
    330344        for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
     
    359373      for (iRow=0;iRow<numberRows;iRow++)
    360374        spare[iRow] = x[iRow]*rowScale[iRow];
    361       if (!hasGaps_) {
     375      if (!(flags_&2)) {
    362376        CoinBigIndex start=columnStart[0];
    363377        for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
     
    419433#endif
    420434  bool packed = rowArray->packedMode();
    421   double factor = 0.27;
     435  double factor = 0.30;
    422436  // We may not want to do by row if there may be cache problems
    423437  // It would be nice to find L2 cache size - for moment 512K
     
    440454  double factor2 = factor*multiplierX;
    441455  if (packed&&rowCopy_&&numberInRowArray>2&&numberInRowArray>factor2*numberRows&&
    442       numberInRowArray<0.9*numberRows) {
     456      numberInRowArray<0.9*numberRows&&0) {
    443457    rowCopy_->transposeTimes(model,rowCopy->getPackedMatrix(),rowArray,y,columnArray);
    444458    return;
     
    447461    // do by column
    448462    // If no gaps - can do a bit faster
    449     if (!hasGaps_) {
     463    if (!(flags_&2)||columnCopy_) {
    450464      transposeTimesByColumn( model,  scalar,
    451465                              rowArray, y, columnArray);
     
    668682        }
    669683      }
     684#if 0
    670685      double value = 0.0;
    671686      CoinBigIndex j;
     
    692707        index[numberNonZero++]=iColumn;
    693708      }
     709#else
     710      if (!columnCopy_)
     711        gutsOfTransposeTimesUnscaled(pi,columnArray,zeroTolerance);
     712      else
     713        columnCopy_->transposeTimes(model,pi,columnArray);
     714      numberNonZero = columnArray->getNumElements();
     715#endif
    694716    } else {
    695717      // scaled
     
    707729      }
    708730      const double * columnScale = model->columnScale();
     731#if 0
    709732      double value = 0.0;
    710733      double scale=columnScale[0];
     
    725748        }
    726749        value = 0.0;
    727         for (j=start; j<end;j++) {
    728           int iRow = row[j];
    729           value += pi[iRow]*elementByColumn[j];
     750        if (model->getColumnStatus(iColumn+1)!=ClpSimplex::basic) {
     751          for (j=start; j<end;j++) {
     752            int iRow = row[j];
     753            value += pi[iRow]*elementByColumn[j];
     754          }
    730755        }
    731756      }
     
    735760        index[numberNonZero++]=iColumn;
    736761      }
     762#else
     763      if (!columnCopy_)
     764        gutsOfTransposeTimesScaled(pi,columnScale,columnArray,zeroTolerance);
     765      else
     766        columnCopy_->transposeTimes(model,pi,columnArray);
     767      numberNonZero = columnArray->getNumElements();
     768#endif
    737769    }
    738770    // zero out
     
    899931    int numberOriginal = 0;
    900932    if (packed) {
     933#if 0
    901934      numberNonZero=0;
    902935      // and set up mark as char array
     
    937970        }
    938971      }
     972#else
     973      gutsOfTransposeTimesByRowGE3(rowArray,columnArray,y,zeroTolerance,scalar);
     974      numberNonZero = columnArray->getNumElements();
     975#endif
    939976    } else {
    940977      double * markVector = y->denseVector();
     
    9821019    double value;
    9831020    if (packed) {
     1021#if 0
    9841022      int iRow0 = whichRow[0];
    9851023      int iRow1 = whichRow[1];
     
    10431081        }
    10441082      }
     1083#else
     1084      gutsOfTransposeTimesByRowEQ2(rowArray,columnArray,y,zeroTolerance,scalar);
     1085      numberNonZero = columnArray->getNumElements();
     1086#endif
    10451087    } else {
    10461088      int iRow = whichRow[0];
     
    10821124    CoinBigIndex j;
    10831125    if (packed) {
     1126#if 0
    10841127      double value = pi[0]*scalar;
    10851128      for (j=rowStart[iRow];j<rowStart[iRow+1];j++) {
     
    10911134        }
    10921135      }
     1136#else
     1137      gutsOfTransposeTimesByRowEQ1(rowArray,columnArray,zeroTolerance,scalar);
     1138      numberNonZero = columnArray->getNumElements();
     1139#endif
    10931140    } else {
    10941141      double value = pi[iRow]*scalar;
     
    11051152  columnArray->setNumElements(numberNonZero);
    11061153  y->setNumElements(0);
     1154}
     1155// Meat of transposeTimes by column when not scaled
     1156void
     1157ClpPackedMatrix::gutsOfTransposeTimesUnscaled(const double * pi,CoinIndexedVector * output, const double zeroTolerance) const
     1158{
     1159  int numberNonZero=0;
     1160  int * index = output->getIndices();
     1161  double * array = output->denseVector();
     1162  // get matrix data pointers
     1163  const int * row = matrix_->getIndices();
     1164  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
     1165  const double * elementByColumn = matrix_->getElements();
     1166  double value = 0.0;
     1167  CoinBigIndex j;
     1168  CoinBigIndex end = columnStart[1];
     1169  for (j=columnStart[0]; j<end;j++) {
     1170    int iRow = row[j];
     1171    value += pi[iRow]*elementByColumn[j];
     1172  }
     1173  int iColumn;
     1174  for (iColumn=0;iColumn<numberActiveColumns_-1;iColumn++) {
     1175    CoinBigIndex start = end;
     1176    end = columnStart[iColumn+2];
     1177    if (fabs(value)>zeroTolerance) {
     1178      array[numberNonZero]=value;
     1179      index[numberNonZero++]=iColumn;
     1180    }
     1181    value = 0.0;
     1182    for (j=start; j<end;j++) {
     1183      int iRow = row[j];
     1184      value += pi[iRow]*elementByColumn[j];
     1185    }
     1186  }
     1187  if (fabs(value)>zeroTolerance) {
     1188    array[numberNonZero]=value;
     1189    index[numberNonZero++]=iColumn;
     1190  }
     1191  output->setNumElements(numberNonZero);
     1192}
     1193// Meat of transposeTimes by column when scaled
     1194void
     1195ClpPackedMatrix::gutsOfTransposeTimesScaled(const double * pi,const double * columnScale,
     1196                                           CoinIndexedVector * output, const double zeroTolerance) const
     1197{
     1198  int numberNonZero=0;
     1199  int * index = output->getIndices();
     1200  double * array = output->denseVector();
     1201  // get matrix data pointers
     1202  const int * row = matrix_->getIndices();
     1203  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
     1204  const double * elementByColumn = matrix_->getElements();
     1205  double value = 0.0;
     1206  double scale=columnScale[0];
     1207  CoinBigIndex j;
     1208  CoinBigIndex end = columnStart[1];
     1209  for (j=columnStart[0]; j<end;j++) {
     1210    int iRow = row[j];
     1211    value += pi[iRow]*elementByColumn[j];
     1212  }
     1213  int iColumn;
     1214  for (iColumn=0;iColumn<numberActiveColumns_-1;iColumn++) {
     1215    value *= scale;
     1216    CoinBigIndex start = end;
     1217    scale = columnScale[iColumn+1];
     1218    end = columnStart[iColumn+2];
     1219    if (fabs(value)>zeroTolerance) {
     1220      array[numberNonZero]=value;
     1221      index[numberNonZero++]=iColumn;
     1222    }
     1223    value = 0.0;
     1224    for (j=start; j<end;j++) {
     1225      int iRow = row[j];
     1226      value += pi[iRow]*elementByColumn[j];
     1227    }
     1228  }
     1229  value *= scale;
     1230  if (fabs(value)>zeroTolerance) {
     1231    array[numberNonZero]=value;
     1232    index[numberNonZero++]=iColumn;
     1233  }
     1234  output->setNumElements(numberNonZero);
     1235}
     1236// Meat of transposeTimes by row n > 2 if packed
     1237void
     1238ClpPackedMatrix::gutsOfTransposeTimesByRowGE3(const CoinIndexedVector * piVector, CoinIndexedVector * output,
     1239                                             CoinIndexedVector * spareVector, const double tolerance, const double scalar) const
     1240{
     1241  double * pi = piVector->denseVector();
     1242  int numberNonZero=0;
     1243  int * index = output->getIndices();
     1244  double * array = output->denseVector();
     1245  int numberInRowArray = piVector->getNumElements();
     1246  const int * column = getIndices();
     1247  const CoinBigIndex * rowStart = getVectorStarts();
     1248  const double * element = getElements();
     1249  const int * whichRow = piVector->getIndices();
     1250  // ** Row copy is already scaled
     1251  int iRow;
     1252  int i;
     1253  // and set up mark as char array
     1254  char * marked = (char *) (index+output->capacity());
     1255  int * lookup = spareVector->getIndices();
     1256#if 1
     1257  for (i=0;i<numberInRowArray;i++) {
     1258    iRow = whichRow[i];
     1259    double value = pi[i]*scalar;
     1260    CoinBigIndex j;
     1261    for (j=rowStart[iRow];j<rowStart[iRow+1];j++) {
     1262      int iColumn = column[j];
     1263      double elValue = element[j];
     1264      if (!marked[iColumn]) {
     1265        marked[iColumn]=1;
     1266        lookup[iColumn]=numberNonZero;
     1267        array[numberNonZero] = value*elValue;
     1268        index[numberNonZero++]=iColumn;
     1269      } else {
     1270        int k = lookup[iColumn];
     1271        array[k] += value*elValue;
     1272      }
     1273    }
     1274  }
     1275#else
     1276  int nextRow = whichRow[0];
     1277  CoinBigIndex nextStart = rowStart[nextRow];
     1278  CoinBigIndex nextEnd = rowStart[nextRow+1];
     1279  whichRow[numberInRowArray]=0; // for electricfence etc
     1280  for (i=0;i<numberInRowArray;i++) {
     1281    iRow = nextRow;
     1282    nextRow = whichRow[i+1];
     1283    CoinBigIndex start = nextStart;
     1284    CoinBigIndex end = nextEnd;
     1285    double value = pi[i]*scalar;
     1286    CoinBigIndex j;
     1287    nextRow = whichRow[i+1];
     1288    nextStart = rowStart[nextRow];
     1289    nextEnd = rowStart[nextRow+1];
     1290    for (j=start;j<end;j++) {
     1291      int iColumn = column[j];
     1292      double elValue = element[j];
     1293      if (!marked[iColumn]) {
     1294        marked[iColumn]=1;
     1295        lookup[iColumn]=numberNonZero;
     1296        array[numberNonZero] = value*elValue;
     1297        index[numberNonZero++]=iColumn;
     1298      } else {
     1299        int k = lookup[iColumn];
     1300        array[k] += value*elValue;
     1301      }
     1302    }
     1303  }
     1304#endif
     1305  // get rid of tiny values and zero out marked
     1306  for (i=0;i<numberNonZero;i++) {
     1307    int iColumn = index[i];
     1308    marked[iColumn]=0;
     1309    double value = array[i];
     1310    while (fabs(value)<=tolerance) {
     1311      numberNonZero--;
     1312      value = array[numberNonZero];
     1313      iColumn = index[numberNonZero];
     1314      marked[iColumn]=0;
     1315      if (i<numberNonZero) {
     1316        array[numberNonZero]=0.0;
     1317        array[i] = value;
     1318        index[i] = iColumn;
     1319      } else {
     1320        array[i]=0.0;
     1321        value =1.0; // to force end of while
     1322      }
     1323    }
     1324  }
     1325  output->setNumElements(numberNonZero);
     1326  spareVector->setNumElements(0);
     1327}
     1328// Meat of transposeTimes by row n == 2 if packed
     1329void
     1330ClpPackedMatrix::gutsOfTransposeTimesByRowEQ2(const CoinIndexedVector * piVector, CoinIndexedVector * output,
     1331                                   CoinIndexedVector * spareVector, const double tolerance, const double scalar) const
     1332{
     1333  double * pi = piVector->denseVector();
     1334  int numberNonZero=0;
     1335  int * index = output->getIndices();
     1336  double * array = output->denseVector();
     1337  const int * column = getIndices();
     1338  const CoinBigIndex * rowStart = getVectorStarts();
     1339  const double * element = getElements();
     1340  const int * whichRow = piVector->getIndices();
     1341  int iRow0 = whichRow[0];
     1342  int iRow1 = whichRow[1];
     1343  double pi0 = pi[0];
     1344  double pi1 = pi[1];
     1345  if (rowStart[iRow0+1]-rowStart[iRow0]>
     1346      rowStart[iRow1+1]-rowStart[iRow1]) {
     1347    // do one with fewer first
     1348    iRow0=iRow1;
     1349    iRow1=whichRow[0];
     1350    pi0=pi1;
     1351    pi1=pi[0];
     1352  }
     1353  // and set up mark as char array
     1354  char * marked = (char *) (index+output->capacity());
     1355  int * lookup = spareVector->getIndices();
     1356  double value = pi0*scalar;
     1357  CoinBigIndex j;
     1358  for (j=rowStart[iRow0];j<rowStart[iRow0+1];j++) {
     1359    int iColumn = column[j];
     1360    double elValue = element[j];
     1361    double value2 = value*elValue;
     1362    array[numberNonZero] = value2;
     1363    marked[iColumn]=1;
     1364    lookup[iColumn]=numberNonZero;
     1365    index[numberNonZero++]=iColumn;
     1366  }
     1367  int numberOriginal = numberNonZero;
     1368  value = pi1*scalar;
     1369  for (j=rowStart[iRow1];j<rowStart[iRow1+1];j++) {
     1370    int iColumn = column[j];
     1371    double elValue = element[j];
     1372    double value2 = value*elValue;
     1373    // I am assuming no zeros in matrix
     1374    if (marked[iColumn]) {
     1375      int iLookup = lookup[iColumn];
     1376      array[iLookup] += value2;
     1377    } else {
     1378      if (fabs(value2)>tolerance) {
     1379        array[numberNonZero] = value2;
     1380        index[numberNonZero++]=iColumn;
     1381      }
     1382    }
     1383  }
     1384  // get rid of tiny values and zero out marked
     1385  int i;
     1386  int iFirst=numberNonZero;
     1387  for (i=0;i<numberOriginal;i++) {
     1388    int iColumn = index[i];
     1389    marked[iColumn]=0;
     1390    if (fabs(array[i])<=tolerance) {
     1391      if (numberNonZero>numberOriginal) {
     1392        numberNonZero--;
     1393        double value = array[numberNonZero];
     1394        array[numberNonZero]=0.0;
     1395        array[i]=value;
     1396        index[i]=index[numberNonZero];
     1397      } else {
     1398        iFirst = i;
     1399      }
     1400    }
     1401  }
     1402 
     1403  if (iFirst<numberNonZero) {
     1404    int n=iFirst;
     1405    for (i=n;i<numberOriginal;i++) {
     1406      int iColumn = index[i];
     1407      double value = array[i];
     1408      array[i]=0.0;
     1409      if (fabs(value)>tolerance) {
     1410        array[n]=value;
     1411        index[n++]=iColumn;
     1412      }
     1413    }
     1414    for (;i<numberNonZero;i++) {
     1415      int iColumn = index[i];
     1416      double value = array[i];
     1417      array[i]=0.0;
     1418      array[n]=value;
     1419      index[n++]=iColumn;
     1420    }
     1421    numberNonZero=n;
     1422  }
     1423  output->setNumElements(numberNonZero);
     1424  spareVector->setNumElements(0);
     1425}
     1426// Meat of transposeTimes by row n == 1 if packed
     1427void
     1428ClpPackedMatrix::gutsOfTransposeTimesByRowEQ1(const CoinIndexedVector * piVector, CoinIndexedVector * output,
     1429                                   const double tolerance, const double scalar) const
     1430{
     1431  double * pi = piVector->denseVector();
     1432  int numberNonZero=0;
     1433  int * index = output->getIndices();
     1434  double * array = output->denseVector();
     1435  const int * column = getIndices();
     1436  const CoinBigIndex * rowStart = getVectorStarts();
     1437  const double * element = getElements();
     1438  int iRow=piVector->getIndices()[0];
     1439  numberNonZero=0;
     1440  CoinBigIndex j;
     1441  double value = pi[0]*scalar;
     1442  for (j=rowStart[iRow];j<rowStart[iRow+1];j++) {
     1443    int iColumn = column[j];
     1444    double elValue = element[j];
     1445    double value2 = value*elValue;
     1446    if (fabs(value2)>tolerance) {
     1447      array[numberNonZero] = value2;
     1448      index[numberNonZero++]=iColumn;
     1449    }
     1450  }
     1451  output->setNumElements(numberNonZero);
    11071452}
    11081453/* Return <code>x *A in <code>z</code> but
     
    11291474  assert (!rowArray->packedMode());
    11301475  columnArray->setPacked();
    1131   if (!hasGaps_&&numberToDo>5) {
     1476  if (!(flags_&2)&&numberToDo>5) {
    11321477    // no gaps and a reasonable number
    11331478    if (!rowScale) {
     
    12221567  bool packed = pi->packedMode();
    12231568  // factor should be smaller if doing both with two pi vectors
    1224   double factor = 0.27;
     1569  double factor = 0.30;
    12251570  // We may not want to do by row if there may be cache problems
    12261571  // It would be nice to find L2 cache size - for moment 512K
     
    12391584  if (!packed)
    12401585    factor *= 0.9;
    1241   return ((numberInRowArray>factor*numberRows||!model->rowCopy())&&!hasGaps_);
     1586  return ((numberInRowArray>factor*numberRows||!model->rowCopy())&&!(flags_&2));
    12421587}
    12431588#ifndef CLP_ALL_ONE_FILE
     
    12911636        pi[iRow]=piOld[i];
    12921637      }
    1293       CoinBigIndex j;
    1294       CoinBigIndex end = columnStart[0];
    1295       for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
    1296         CoinBigIndex start = end;
    1297         end = columnStart[iColumn+1];
    1298         ClpSimplex::Status status = model->getStatus(iColumn);
    1299         if (status==ClpSimplex::basic||status==ClpSimplex::isFixed) continue;
    1300         double value = 0.0;
    1301         for (j=start; j<end;j++) {
    1302           int iRow = row[j];
    1303           value -= pi[iRow]*elementByColumn[j];
    1304         }
    1305         if (fabs(value)>zeroTolerance) {
    1306           // and do other array
    1307           double modification = 0.0;
    1308           for (j=start; j<end;j++) {
    1309             int iRow = row[j];
    1310             modification += piWeight[iRow]*elementByColumn[j];
    1311           }
    1312           double thisWeight = weights[iColumn];
    1313           double pivot = value*scaleFactor;
    1314           double pivotSquared = pivot * pivot;
    1315           thisWeight += pivotSquared * devex + pivot * modification;
    1316           if (thisWeight<DEVEX_TRY_NORM) {
    1317             if (referenceIn<0.0) {
    1318               // steepest
    1319               thisWeight = CoinMax(DEVEX_TRY_NORM,DEVEX_ADD_ONE+pivotSquared);
    1320             } else {
    1321               // exact
    1322               thisWeight = referenceIn*pivotSquared;
    1323               if (reference(iColumn))
    1324                 thisWeight += 1.0;
    1325               thisWeight = CoinMax(thisWeight,DEVEX_TRY_NORM);
    1326             }
    1327           }
    1328           weights[iColumn] = thisWeight;
    1329           if (!killDjs) {
    1330             array[numberNonZero]=value;
    1331             index[numberNonZero++]=iColumn;
    1332           }
    1333         }
     1638      if (!columnCopy_) {
     1639        CoinBigIndex j;
     1640        CoinBigIndex end = columnStart[0];
     1641        for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
     1642          CoinBigIndex start = end;
     1643          end = columnStart[iColumn+1];
     1644          ClpSimplex::Status status = model->getStatus(iColumn);
     1645          if (status==ClpSimplex::basic||status==ClpSimplex::isFixed) continue;
     1646          double value = 0.0;
     1647          for (j=start; j<end;j++) {
     1648            int iRow = row[j];
     1649            value -= pi[iRow]*elementByColumn[j];
     1650          }
     1651          if (fabs(value)>zeroTolerance) {
     1652            // and do other array
     1653            double modification = 0.0;
     1654            for (j=start; j<end;j++) {
     1655              int iRow = row[j];
     1656              modification += piWeight[iRow]*elementByColumn[j];
     1657            }
     1658            double thisWeight = weights[iColumn];
     1659            double pivot = value*scaleFactor;
     1660            double pivotSquared = pivot * pivot;
     1661            thisWeight += pivotSquared * devex + pivot * modification;
     1662            if (thisWeight<DEVEX_TRY_NORM) {
     1663              if (referenceIn<0.0) {
     1664                // steepest
     1665                thisWeight = CoinMax(DEVEX_TRY_NORM,DEVEX_ADD_ONE+pivotSquared);
     1666              } else {
     1667                // exact
     1668                thisWeight = referenceIn*pivotSquared;
     1669                if (reference(iColumn))
     1670                  thisWeight += 1.0;
     1671                thisWeight = CoinMax(thisWeight,DEVEX_TRY_NORM);
     1672              }
     1673            }
     1674            weights[iColumn] = thisWeight;
     1675            if (!killDjs) {
     1676              array[numberNonZero]=value;
     1677              index[numberNonZero++]=iColumn;
     1678            }
     1679          }
     1680        }
     1681      } else {
     1682        // use special column copy
     1683        // reset back
     1684        if (killDjs)
     1685          scaleFactor=0.0;
     1686        columnCopy_->transposeTimes2(model,pi,dj1,piWeight,referenceIn,devex,
     1687                                     reference,weights,scaleFactor);
     1688        numberNonZero = dj1->getNumElements();
    13341689      }
    13351690    } else {
     
    13401695        pi[iRow]=piOld[i]*rowScale[iRow];
    13411696      }
    1342       const double * columnScale = model->columnScale();
    1343       CoinBigIndex j;
    1344       CoinBigIndex end = columnStart[0];
    1345       for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
    1346         CoinBigIndex start = end;
    1347         end = columnStart[iColumn+1];
    1348         ClpSimplex::Status status = model->getStatus(iColumn);
    1349         if (status==ClpSimplex::basic||status==ClpSimplex::isFixed) continue;
    1350         double scale=columnScale[iColumn];
    1351         double value = 0.0;
    1352         for (j=start; j<end;j++) {
    1353           int iRow = row[j];
    1354           value -= pi[iRow]*elementByColumn[j];
    1355         }
    1356         value *= scale;
    1357         if (fabs(value)>zeroTolerance) {
    1358           double modification = 0.0;
    1359           for (j=start; j<end;j++) {
    1360             int iRow = row[j];
    1361             modification += piWeight[iRow]*elementByColumn[j]*rowScale[iRow];
    1362           }
    1363           modification *= scale;
    1364           double thisWeight = weights[iColumn];
    1365           double pivot = value*scaleFactor;
    1366           double pivotSquared = pivot * pivot;
    1367           thisWeight += pivotSquared * devex + pivot * modification;
    1368           if (thisWeight<DEVEX_TRY_NORM) {
    1369             if (referenceIn<0.0) {
    1370               // steepest
    1371               thisWeight = CoinMax(DEVEX_TRY_NORM,DEVEX_ADD_ONE+pivotSquared);
    1372             } else {
    1373               // exact
    1374               thisWeight = referenceIn*pivotSquared;
    1375               if (reference(iColumn))
    1376                 thisWeight += 1.0;
    1377               thisWeight = CoinMax(thisWeight,DEVEX_TRY_NORM);
    1378             }
    1379           }
    1380           weights[iColumn] = thisWeight;
    1381           if (!killDjs) {
    1382             array[numberNonZero]=value;
    1383             index[numberNonZero++]=iColumn;
    1384           }
    1385         }
     1697      // can also scale piWeight as not used again
     1698      int numberWeight = pi2->getNumElements();
     1699      const int * indexWeight = pi2->getIndices();
     1700      for (i=0;i<numberWeight;i++) {
     1701        int iRow = indexWeight[i];
     1702        piWeight[iRow] *= rowScale[iRow];
     1703      }
     1704      if (!columnCopy_) {
     1705        const double * columnScale = model->columnScale();
     1706        CoinBigIndex j;
     1707        CoinBigIndex end = columnStart[0];
     1708        for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
     1709          CoinBigIndex start = end;
     1710          end = columnStart[iColumn+1];
     1711          ClpSimplex::Status status = model->getStatus(iColumn);
     1712          if (status==ClpSimplex::basic||status==ClpSimplex::isFixed) continue;
     1713          double scale=columnScale[iColumn];
     1714          double value = 0.0;
     1715          for (j=start; j<end;j++) {
     1716            int iRow = row[j];
     1717            value -= pi[iRow]*elementByColumn[j];
     1718          }
     1719          value *= scale;
     1720          if (fabs(value)>zeroTolerance) {
     1721            double modification = 0.0;
     1722            for (j=start; j<end;j++) {
     1723              int iRow = row[j];
     1724              modification += piWeight[iRow]*elementByColumn[j];
     1725            }
     1726            modification *= scale;
     1727            double thisWeight = weights[iColumn];
     1728            double pivot = value*scaleFactor;
     1729            double pivotSquared = pivot * pivot;
     1730            thisWeight += pivotSquared * devex + pivot * modification;
     1731            if (thisWeight<DEVEX_TRY_NORM) {
     1732              if (referenceIn<0.0) {
     1733                // steepest
     1734                thisWeight = CoinMax(DEVEX_TRY_NORM,DEVEX_ADD_ONE+pivotSquared);
     1735              } else {
     1736                // exact
     1737                thisWeight = referenceIn*pivotSquared;
     1738                if (reference(iColumn))
     1739                  thisWeight += 1.0;
     1740                thisWeight = CoinMax(thisWeight,DEVEX_TRY_NORM);
     1741              }
     1742            }
     1743            weights[iColumn] = thisWeight;
     1744            if (!killDjs) {
     1745              array[numberNonZero]=value;
     1746              index[numberNonZero++]=iColumn;
     1747            }
     1748          }
     1749        }
     1750      } else {
     1751        // use special column copy
     1752        // reset back
     1753        if (killDjs)
     1754          scaleFactor=0.0;
     1755        columnCopy_->transposeTimes2(model,pi,dj1,piWeight,referenceIn,devex,
     1756                                     reference,weights,scaleFactor);
     1757        numberNonZero = dj1->getNumElements();
    13861758      }
    13871759    }
     
    14421814    } else {
    14431815      // scaled
     1816      // can also scale piWeight as not used again
     1817      int numberWeight = pi2->getNumElements();
     1818      const int * indexWeight = pi2->getIndices();
     1819      for (int i=0;i<numberWeight;i++) {
     1820        int iRow = indexWeight[i];
     1821        piWeight[iRow] *= rowScale[iRow];
     1822      }
    14441823      const double * columnScale = model->columnScale();
    14451824      CoinBigIndex j;
     
    14611840          for (j=start; j<end;j++) {
    14621841            int iRow = row[j];
    1463             modification += piWeight[iRow]*elementByColumn[j]*rowScale[iRow];
     1842            modification += piWeight[iRow]*elementByColumn[j];
    14641843          }
    14651844          modification *= scale;
     
    16181997  const int * row = matrix_->getIndices();
    16191998  const double * elementByColumn = matrix_->getElements();
    1620   if (!zeroElements_) {
     1999  if ((flags_&1)==0) {
    16212000    if (!rowScale) {
    16222001      // no scaling
     
    17002079ClpPackedMatrix::scale(ClpModel * model) const
    17012080{
     2081#ifndef NDEBUG
     2082  //checkFlags();
     2083#endif
    17022084  int numberRows = model->numberRows();
    17032085  int numberColumns = matrix_->getNumCols();
     
    21772559  const double * elementByColumn = matrix_->getElements();
    21782560  CoinBigIndex i;
     2561  int * index = rowArray->getIndices();
     2562  double * array = rowArray->denseVector();
     2563  int number = 0;
    21792564  if (!rowScale) {
    2180     CoinBigIndex j=columnStart[iColumn];
    2181     rowArray->createPacked(columnLength[iColumn],
    2182                            row+j,elementByColumn+j);
     2565    for (i=columnStart[iColumn];
     2566         i<columnStart[iColumn]+columnLength[iColumn];i++) {
     2567      int iRow = row[i];
     2568      double value = elementByColumn[i];
     2569      if (value) {
     2570        array[number]=value;
     2571        index[number++]=iRow;
     2572      }
     2573    }
     2574    rowArray->setNumElements(number);
     2575    rowArray->setPackedMode(true);
    21832576  } else {
    21842577    // apply scaling
    21852578    double scale = model->columnScale()[iColumn];
    2186     int * index = rowArray->getIndices();
    2187     double * array = rowArray->denseVector();
    2188     int number = 0;
    21892579    for (i=columnStart[iColumn];
    21902580         i<columnStart[iColumn]+columnLength[iColumn];i++) {
    21912581      int iRow = row[i];
    2192       array[number]=elementByColumn[i]*scale*rowScale[iRow];
    2193       index[number++]=iRow;
     2582      double value = elementByColumn[i]*scale*rowScale[iRow];
     2583      if (value) {
     2584        array[number]=value;
     2585        index[number++]=iRow;
     2586      }
    21942587    }
    21952588    rowArray->setNumElements(number);
     
    22842677  int numberColumns = matrix_->getNumCols();
    22852678  // Say no gaps
    2286   hasGaps_=false;
     2679  flags_ &= ~2;
    22872680  if (check==14) {
    22882681    for (iColumn=0;iColumn<numberColumns;iColumn++) {
     
    22902683      CoinBigIndex end = start + columnLength[iColumn];
    22912684      if (end!=columnStart[iColumn+1]) {
    2292         hasGaps_=true;
     2685        flags_ |= 2;
    22932686        break;
    22942687      }
     
    23062699    CoinBigIndex end = start + columnLength[iColumn];
    23072700    if (end!=columnStart[iColumn+1])
    2308       hasGaps_=true;
     2701      flags_ |= 2;
    23092702    for (j=start;j<end;j++) {
    23102703      double value = fabs(elementByColumn[j]);
     
    23282721      //printf("%d %d %d %g\n",iColumn,j,row[j],elementByColumn[j]);
    23292722      if (!value)
    2330         zeroElements_ = true; // there are zero elements
     2723        flags_ |= 1; // there are zero elements
    23312724      if (value<smallest) {
    23322725        numberSmall++;
     
    23692762  // If smallest >0.0 then there can't be zero elements
    23702763  if (smallest>0.0)
    2371     zeroElements_=false;
     2764    flags_ &= ~1;;
    23722765  return true;
    23732766}
     
    26883081  rhsOffset(model,true);
    26893082}
     3083// Gets rid of special copies
     3084void
     3085ClpPackedMatrix::clearCopies()
     3086{
     3087  delete rowCopy_;
     3088  delete columnCopy_;
     3089  rowCopy_=NULL;
     3090  columnCopy_=NULL;
     3091  flags_ &= ~(4+8);
     3092 }
    26903093// makes sure active columns correct
    26913094int
     
    26933096{
    26943097  numberActiveColumns_ = matrix_->getNumCols();
     3098#if 0
     3099  ClpMatrixBase * rowCopyBase = reverseOrderedCopy();
     3100  ClpPackedMatrix* rowCopy =
     3101    dynamic_cast< ClpPackedMatrix*>(rowCopyBase);
     3102  // Make sure it is really a ClpPackedMatrix
     3103  assert (rowCopy!=NULL);
     3104
     3105  const int * column = rowCopy->getIndices();
     3106  const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
     3107  const int * rowLength = rowCopy->getVectorLengths();
     3108  const double * element = rowCopy->getElements();
     3109  for (int i=0;i<rowCopy->getNumRows();i++) {
     3110    if (!rowLength[i])
     3111      printf("zero row %d\n",i);
     3112  }
     3113  delete rowCopy;
     3114#endif
    26953115  return 0;
    26963116}
     
    27773197ClpPackedMatrix::reallyScale(const double * rowScale, const double * columnScale)
    27783198{
     3199  clearCopies();
    27793200  int numberColumns = matrix_->getNumCols();
    27803201  const int * row = matrix_->getIndices();
     
    27983219  if (matrix_->getNumCols())
    27993220    matrix_->deleteCols(numDel,indDel);
     3221  clearCopies();
    28003222  numberActiveColumns_ = matrix_->getNumCols();
    28013223  // may now have gaps
    2802   hasGaps_=true;
     3224  flags_ |= 2;
    28033225  matrix_->setExtraGap(1.0e-50);
    28043226}
     
    28093231  if (matrix_->getNumRows())
    28103232    matrix_->deleteRows(numDel,indDel);
     3233  clearCopies();
    28113234  numberActiveColumns_ = matrix_->getNumCols();
    28123235  // may now have gaps
    2813   hasGaps_=true;
     3236  flags_ |= 2;
    28143237  matrix_->setExtraGap(1.0e-50);
    28153238}
     
    28213244  matrix_->appendCols(number,columns);
    28223245  numberActiveColumns_ = matrix_->getNumCols();
     3246  clearCopies();
    28233247}
    28243248// Append Rows
     
    28293253  numberActiveColumns_ = matrix_->getNumCols();
    28303254  // may now have gaps
    2831   hasGaps_=true;
     3255  flags_ |= 2;
     3256  clearCopies();
    28323257}
    28333258#endif
     
    28633288      matrix_->setDimensions(numberOther,-1);
    28643289    numberErrors=matrix_->appendCols(number,starts,index,element,numberOther);
    2865     numberActiveColumns_ = matrix_->getNumCols();
    2866   }
     3290  }
     3291  clearCopies();
    28673292  return numberErrors;
    28683293}
     
    28763301    delete rowCopy_;
    28773302    rowCopy_=NULL;
     3303    flags_ &= ~4;
     3304  } else {
     3305    flags_ |= 4;
     3306  }
     3307}
     3308void
     3309ClpPackedMatrix::specialColumnCopy(ClpSimplex * model)
     3310{
     3311  delete columnCopy_;
     3312  if ((flags_&8)!=0)
     3313    columnCopy_ = new ClpPackedMatrix3(model,matrix_);
     3314  else
     3315    columnCopy_=NULL;
     3316}
     3317// Say we don't want special column copy
     3318void
     3319ClpPackedMatrix::releaseSpecialColumnCopy()
     3320{
     3321  flags_ &= ~8;
     3322  delete columnCopy_;
     3323  columnCopy_=NULL;
     3324}
     3325// Correct sequence in and out to give true value
     3326void
     3327ClpPackedMatrix::correctSequence(const ClpSimplex * model,int & sequenceIn, int & sequenceOut)
     3328{
     3329  if (columnCopy_) {
     3330    if (sequenceIn!=-999) {
     3331      if (sequenceIn!=sequenceOut) {
     3332        if (sequenceIn<numberActiveColumns_)
     3333          columnCopy_->swapOne(model,this,sequenceIn);
     3334        if (sequenceOut<numberActiveColumns_)
     3335          columnCopy_->swapOne(model,this,sequenceOut);
     3336      }
     3337    } else {
     3338      // do all
     3339      columnCopy_->sortBlocks(model);
     3340    }
     3341  }
     3342}
     3343// Check validity
     3344void
     3345ClpPackedMatrix::checkFlags() const
     3346{
     3347  int iColumn;
     3348  // get matrix data pointers
     3349  //const int * row = matrix_->getIndices();
     3350  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
     3351  const int * columnLength = matrix_->getVectorLengths();
     3352  const double * elementByColumn = matrix_->getElements();
     3353  if (!zeros()) {
     3354    for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
     3355      CoinBigIndex j;
     3356      for (j=columnStart[iColumn];j<columnStart[iColumn]+columnLength[iColumn];j++) {
     3357        if (!elementByColumn[j])
     3358          abort();
     3359      }
     3360    }
     3361  }
     3362  if ((flags_&2)==0) {
     3363    for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
     3364      if (columnStart[iColumn+1]!=columnStart[iColumn]+columnLength[iColumn]) {
     3365        abort();
     3366      }
     3367    }
    28783368  }
    28793369}
     
    29273417  // tune
    29283418#if 0
    2929   int chunkY[5]={4096,8192,16384,32768,65535};
     3419  int chunkY[7]={1024,2048,4096,8192,16384,32768,65535};
    29303420  int its = model->maximumIterations();
    29313421  if (its>=1000000&&its<1000999) {
    29323422    its -= 1000000;
    29333423    its = its/10;
    2934     if (its>=5) abort();
     3424    if (its>=7) abort();
    29353425    chunk=chunkY[its];
    29363426    printf("chunk size %d\n",chunk);
     
    37174207  }
    37184208}
     4209/* Default constructor. */
     4210ClpPackedMatrix3::ClpPackedMatrix3()
     4211  : numberBlocks_(0),
     4212    numberColumns_(0),
     4213    column_(NULL),
     4214    start_(NULL),
     4215    row_(NULL),
     4216    element_(NULL),
     4217    block_(NULL)
     4218{
     4219}
     4220/* Constructor from copy. */
     4221ClpPackedMatrix3::ClpPackedMatrix3(ClpSimplex * model,const CoinPackedMatrix * columnCopy)
     4222  : numberBlocks_(0),
     4223    numberColumns_(0),
     4224    column_(NULL),
     4225    start_(NULL),
     4226    row_(NULL),
     4227    element_(NULL),
     4228    block_(NULL)
     4229{
     4230#define MINBLOCK 6
     4231#define MAXBLOCK 100
     4232#define MAXUNROLL 10
     4233  numberColumns_ = columnCopy->getNumCols();
     4234  int numberRows = columnCopy->getNumRows();
     4235  int * counts = new int[numberRows+1];
     4236  CoinZeroN(counts,numberRows+1);
     4237  CoinBigIndex nels=0;
     4238  CoinBigIndex nZeroEl=0;
     4239  int iColumn;
     4240  // get matrix data pointers
     4241  const int * row = columnCopy->getIndices();
     4242  const CoinBigIndex * columnStart = columnCopy->getVectorStarts();
     4243  const int * columnLength = columnCopy->getVectorLengths();
     4244  const double * elementByColumn = columnCopy->getElements();
     4245  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
     4246    CoinBigIndex start = columnStart[iColumn];
     4247    int n = columnLength[iColumn];
     4248    CoinBigIndex end = start + n;
     4249    int kZero=0;
     4250    for (CoinBigIndex j=start;j<end;j++) {
     4251      if(!elementByColumn[j])
     4252        kZero++;
     4253    }
     4254    n -= kZero;
     4255    nZeroEl += kZero;
     4256    nels += n;
     4257    counts[n]++;
     4258  }
     4259  int nZeroColumns = counts[0];
     4260  counts[0]=-1;
     4261  numberColumns_ -= nZeroColumns;
     4262  column_ = new int [2*numberColumns_];
     4263  int * lookup = column_ + numberColumns_;
     4264  row_ = new int[nels];
     4265  element_ = new double[nels];
     4266  int nOdd=0;
     4267  CoinBigIndex nInOdd=0;
     4268  int i;
     4269  for (i=1;i<=numberRows;i++) {
     4270    int n = counts[i];
     4271    if (n) {
     4272      if (n<MINBLOCK||i>MAXBLOCK) {
     4273        nOdd += n;
     4274        counts[i]=-1;
     4275        nInOdd += n*i;
     4276      } else {
     4277        numberBlocks_++;
     4278      }
     4279    } else {
     4280      counts[i]=-1;
     4281    }
     4282  }
     4283  start_ = new CoinBigIndex [nOdd+1];
     4284  // even if no blocks do a dummy one
     4285  numberBlocks_ = CoinMax(numberBlocks_,1);
     4286  block_ = (blockStruct *) new char [numberBlocks_*sizeof(blockStruct)];
     4287  memset(block_,0,numberBlocks_*sizeof(blockStruct));
     4288  // Fill in what we can
     4289  int nTotal=nOdd;
     4290  // in case no blocks
     4291  block_->startIndices_=nTotal;
     4292  nels=nInOdd;
     4293  int nBlock=0;
     4294  for (i=0;i<=CoinMin(MAXBLOCK,numberRows);i++) {
     4295    if (counts[i]>0) {
     4296      blockStruct * block = block_ + nBlock;
     4297      int n=counts[i];
     4298      counts[i]= nBlock; // backward pointer
     4299      nBlock++;
     4300      block->startIndices_ = nTotal;
     4301      block->startElements_ = nels;
     4302      block->numberElements_ = i;
     4303      // up counts
     4304      nTotal += n;
     4305      nels += n*i;
     4306    }
     4307  }
     4308  // fill
     4309  start_[0]=0;
     4310  nOdd=0;
     4311  nInOdd=0;
     4312  const double * columnScale = model->columnScale();
     4313  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
     4314    CoinBigIndex start = columnStart[iColumn];
     4315    int n = columnLength[iColumn];
     4316    CoinBigIndex end = start + n;
     4317    int kZero=0;
     4318    for (CoinBigIndex j=start;j<end;j++) {
     4319      if(!elementByColumn[j])
     4320        kZero++;
     4321    }
     4322    n -= kZero;
     4323    if (n) {
     4324      int iBlock = counts[n];
     4325      if (iBlock>=0) {
     4326        blockStruct * block = block_ + iBlock;
     4327        int k = block->numberInBlock_;
     4328        block->numberInBlock_ ++;
     4329        column_[block->startIndices_+k]=iColumn;
     4330        lookup[iColumn]=k;
     4331        CoinBigIndex put = block->startElements_+k*n;
     4332        for (CoinBigIndex j=start;j<end;j++) {
     4333          double value = elementByColumn[j];
     4334          if(value) {
     4335            if (columnScale)
     4336              value *= columnScale[iColumn];
     4337            element_[put] = value;
     4338            row_[put++]=row[j];
     4339          }
     4340        }
     4341      } else {
     4342        // odd ones
     4343        for (CoinBigIndex j=start;j<end;j++) {
     4344          double value = elementByColumn[j];
     4345          if(value) {
     4346            if (columnScale)
     4347              value *= columnScale[iColumn];
     4348            element_[nInOdd] = value;
     4349            row_[nInOdd++]=row[j];
     4350          }
     4351        }
     4352        column_[nOdd]=iColumn;
     4353        lookup[iColumn]=-1;
     4354        nOdd++;
     4355        start_[nOdd]=nInOdd;
     4356      }
     4357    }
     4358  }
     4359  delete [] counts;
     4360}
     4361/* Destructor */
     4362ClpPackedMatrix3::~ClpPackedMatrix3()
     4363{
     4364  delete [] column_;
     4365  delete [] start_;
     4366  delete [] row_;
     4367  delete [] element_;
     4368  delete [] block_;
     4369}
     4370/* The copy constructor. */
     4371ClpPackedMatrix3::ClpPackedMatrix3(const ClpPackedMatrix3 & rhs)
     4372  : numberBlocks_(rhs.numberBlocks_),
     4373    numberColumns_(rhs.numberColumns_),
     4374    column_(NULL),
     4375    start_(NULL),
     4376    row_(NULL),
     4377    element_(NULL),
     4378    block_(NULL)
     4379{
     4380  if (rhs.numberBlocks_) {
     4381    block_ = CoinCopyOfArray(rhs.block_,numberBlocks_);
     4382    column_ = CoinCopyOfArray(rhs.column_,2*numberColumns_);
     4383    int numberOdd = block_->startIndices_;
     4384    start_ = CoinCopyOfArray(rhs.start_,numberOdd+1);
     4385    blockStruct * lastBlock = block_ + (numberBlocks_-1);
     4386    CoinBigIndex numberElements = lastBlock->startElements_+
     4387      lastBlock->numberInBlock_*lastBlock->numberElements_;
     4388    row_ = CoinCopyOfArray(rhs.row_,numberElements);
     4389    element_ = CoinCopyOfArray(rhs.element_,numberElements);
     4390  }
     4391}
     4392ClpPackedMatrix3&
     4393ClpPackedMatrix3::operator=(const ClpPackedMatrix3 & rhs)
     4394{
     4395  if (this != &rhs) {
     4396    delete [] column_;
     4397    delete [] start_;
     4398    delete [] row_;
     4399    delete [] element_;
     4400    delete [] block_;
     4401    numberBlocks_ = rhs.numberBlocks_;
     4402    numberColumns_ = rhs.numberColumns_;
     4403    if (rhs.numberBlocks_) {
     4404      block_ = CoinCopyOfArray(rhs.block_,numberBlocks_);
     4405      column_ = CoinCopyOfArray(rhs.column_,2*numberColumns_);
     4406      int numberOdd = block_->startIndices_;
     4407      start_ = CoinCopyOfArray(rhs.start_,numberOdd+1);
     4408      blockStruct * lastBlock = block_ + (numberBlocks_-1);
     4409      CoinBigIndex numberElements = lastBlock->startElements_+
     4410        lastBlock->numberInBlock_*lastBlock->numberElements_;
     4411      row_ = CoinCopyOfArray(rhs.row_,numberElements);
     4412      element_ = CoinCopyOfArray(rhs.element_,numberElements);
     4413    } else {
     4414      column_ = NULL;
     4415      start_ = NULL;
     4416      row_ = NULL;
     4417      element_ = NULL;
     4418      block_ = NULL;
     4419    }
     4420  }
     4421  return *this;
     4422}
     4423/* Sort blocks */
     4424void
     4425ClpPackedMatrix3::sortBlocks(const ClpSimplex * model)
     4426{
     4427  int * lookup = column_ + numberColumns_;
     4428  for (int iBlock=0;iBlock<numberBlocks_;iBlock++) {
     4429    blockStruct * block = block_ + iBlock;
     4430    int numberInBlock = block->numberInBlock_;
     4431    int nel = block->numberElements_;
     4432    int * row = row_ + block->startElements_;
     4433    double * element = element_ + block->startElements_;
     4434    int * column = column_ + block->startIndices_;
     4435    int lastPrice=0;
     4436    int firstNotPrice=numberInBlock-1;
     4437    while (lastPrice<=firstNotPrice) {
     4438      // find first basic or fixed
     4439      int iColumn=numberInBlock;
     4440      for (;lastPrice<=firstNotPrice;lastPrice++) {
     4441        iColumn = column[lastPrice];
     4442        if (model->getColumnStatus(iColumn)==ClpSimplex::basic||
     4443            model->getColumnStatus(iColumn)==ClpSimplex::isFixed)
     4444          break;
     4445      }
     4446      // find last non basic or fixed
     4447      int jColumn=-1;
     4448      for (;firstNotPrice>lastPrice;firstNotPrice--) {
     4449        jColumn = column[firstNotPrice];
     4450        if (model->getColumnStatus(jColumn)!=ClpSimplex::basic&&
     4451            model->getColumnStatus(jColumn)!=ClpSimplex::isFixed)
     4452          break;
     4453      }
     4454      if (firstNotPrice>lastPrice) {
     4455        assert (column[lastPrice]==iColumn);
     4456        assert (column[firstNotPrice]==jColumn);
     4457        // need to swap
     4458        column[firstNotPrice]=iColumn;
     4459        lookup[iColumn]=firstNotPrice;
     4460        column[lastPrice]=jColumn;
     4461        lookup[jColumn]=lastPrice;
     4462        double * elementA = element + lastPrice*nel;
     4463        int * rowA = row + lastPrice*nel;
     4464        double * elementB = element + firstNotPrice*nel;
     4465        int * rowB = row + firstNotPrice*nel;
     4466        for (int i=0;i<nel;i++) {
     4467          int temp = rowA[i];
     4468          double tempE = elementA[i];
     4469          rowA[i]=rowB[i];
     4470          elementA[i]=elementB[i];
     4471          rowB[i]=temp;
     4472          elementB[i]=tempE;
     4473        }
     4474        firstNotPrice--;
     4475        lastPrice++;
     4476      } else if (lastPrice==firstNotPrice) {
     4477        // make sure correct side
     4478        iColumn = column[lastPrice];
     4479        if (model->getColumnStatus(iColumn)!=ClpSimplex::basic&&
     4480            model->getColumnStatus(iColumn)!=ClpSimplex::isFixed)
     4481          lastPrice++;
     4482        break;
     4483      }
     4484    }
     4485    block->numberPrice_=lastPrice;
     4486#ifndef NDEBUG
     4487    // check
     4488    int i;
     4489    for (i=0;i<lastPrice;i++) {
     4490      int iColumn = column[i];
     4491      assert (model->getColumnStatus(iColumn)!=ClpSimplex::basic&&
     4492              model->getColumnStatus(iColumn)!=ClpSimplex::isFixed);
     4493      assert (lookup[iColumn]==i);
     4494    }
     4495    for (;i<numberInBlock;i++) {
     4496      int iColumn = column[i];
     4497      assert (model->getColumnStatus(iColumn)==ClpSimplex::basic||
     4498              model->getColumnStatus(iColumn)==ClpSimplex::isFixed);
     4499      assert (lookup[iColumn]==i);
     4500    }
     4501#endif
     4502  }
     4503}
     4504// Swap one variable
     4505void
     4506ClpPackedMatrix3::swapOne(const ClpSimplex * model,const ClpPackedMatrix * matrix,
     4507                          int iColumn)
     4508{
     4509  int * lookup = column_ + numberColumns_;
     4510  // position in block
     4511  int kA = lookup[iColumn];
     4512  if (kA<0)
     4513    return; // odd one
     4514  // get matrix data pointers
     4515  const CoinPackedMatrix * columnCopy = matrix->getPackedMatrix();
     4516  //const int * row = columnCopy->getIndices();
     4517  const CoinBigIndex * columnStart = columnCopy->getVectorStarts();
     4518  const int * columnLength = columnCopy->getVectorLengths();
     4519  const double * elementByColumn = columnCopy->getElements();
     4520  CoinBigIndex start = columnStart[iColumn];
     4521  int n = columnLength[iColumn];
     4522  if (matrix->zeros()) {
     4523    CoinBigIndex end = start + n;
     4524    for (CoinBigIndex j=start;j<end;j++) {
     4525      if(!elementByColumn[j])
     4526        n--;
     4527    }
     4528  }
     4529  // find block - could do binary search
     4530  int iBlock=CoinMin(n,numberBlocks_)-1;
     4531  while (block_[iBlock].numberElements_!=n)
     4532    iBlock--;
     4533  blockStruct * block = block_ + iBlock;
     4534  int nel = block->numberElements_;
     4535  int * row = row_ + block->startElements_;
     4536  double * element = element_ + block->startElements_;
     4537  int * column = column_ + block->startIndices_;
     4538  assert (column[kA]==iColumn);
     4539  bool moveUp = (model->getColumnStatus(iColumn)==ClpSimplex::basic||
     4540                 model->getColumnStatus(iColumn)==ClpSimplex::isFixed);
     4541  int lastPrice=block->numberPrice_;
     4542  int kB;
     4543  if (moveUp) {
     4544    kB=lastPrice-1;
     4545    block->numberPrice_--;
     4546  } else {
     4547    kB=lastPrice;
     4548    block->numberPrice_++;
     4549  }
     4550  int jColumn = column[kB];
     4551  column[kA]=jColumn;
     4552  lookup[jColumn]=kA;
     4553  column[kB]=iColumn;
     4554  lookup[iColumn]=kB;
     4555  double * elementA = element + kB*nel;
     4556  int * rowA = row + kB*nel;
     4557  double * elementB = element + kA*nel;
     4558  int * rowB = row + kA*nel;
     4559  for (int i=0;i<nel;i++) {
     4560    int temp = rowA[i];
     4561    double tempE = elementA[i];
     4562    rowA[i]=rowB[i];
     4563    elementA[i]=elementB[i];
     4564    rowB[i]=temp;
     4565    elementB[i]=tempE;
     4566  }
     4567#if 1
     4568#ifndef NDEBUG
     4569  // check
     4570  int i;
     4571  for (i=0;i<block->numberPrice_;i++) {
     4572    int iColumn = column[i];
     4573    if (iColumn!=model->sequenceIn()&&iColumn!=model->sequenceOut())
     4574      assert (model->getColumnStatus(iColumn)!=ClpSimplex::basic&&
     4575              model->getColumnStatus(iColumn)!=ClpSimplex::isFixed);
     4576    assert (lookup[iColumn]==i);
     4577  }
     4578  int numberInBlock = block->numberInBlock_;
     4579  for (;i<numberInBlock;i++) {
     4580    int iColumn = column[i];
     4581    if (iColumn!=model->sequenceIn()&&iColumn!=model->sequenceOut())
     4582      assert (model->getColumnStatus(iColumn)==ClpSimplex::basic||
     4583              model->getColumnStatus(iColumn)==ClpSimplex::isFixed);
     4584    assert (lookup[iColumn]==i);
     4585  }
     4586#endif
     4587#endif
     4588}
     4589/* Return <code>x * -1 * A in <code>z</code>.
     4590   Note - x packed and z will be packed mode
     4591   Squashes small elements and knows about ClpSimplex */
     4592void
     4593ClpPackedMatrix3::transposeTimes(const ClpSimplex * model,
     4594                                 const double * pi,
     4595                                 CoinIndexedVector * output) const
     4596{
     4597  int numberNonZero=0;
     4598  int * index = output->getIndices();
     4599  double * array = output->denseVector();
     4600  double zeroTolerance = model->factorization()->zeroTolerance();
     4601  double value = 0.0;
     4602  CoinBigIndex j;
     4603  int numberOdd = block_->startIndices_;
     4604  if (numberOdd) {
     4605    // A) as probably long may be worth unrolling
     4606    CoinBigIndex end = start_[1];
     4607    for (j=start_[0]; j<end;j++) {
     4608      int iRow = row_[j];
     4609      value += pi[iRow]*element_[j];
     4610    }
     4611    int iColumn;
     4612    // int jColumn=column_[0];
     4613   
     4614    for (iColumn=0;iColumn<numberOdd-1;iColumn++) {
     4615      CoinBigIndex start = end;
     4616      end = start_[iColumn+2];
     4617      if (fabs(value)>zeroTolerance) {
     4618        array[numberNonZero]=value;
     4619        index[numberNonZero++]=column_[iColumn];
     4620        //index[numberNonZero++]=jColumn;
     4621      }
     4622      // jColumn = column_[iColumn+1];
     4623      value = 0.0;
     4624      //if (model->getColumnStatus(jColumn)!=ClpSimplex::basic) {
     4625        for (j=start; j<end;j++) {
     4626          int iRow = row_[j];
     4627          value += pi[iRow]*element_[j];
     4628        }
     4629        //}
     4630    }
     4631    if (fabs(value)>zeroTolerance) {
     4632      array[numberNonZero]=value;
     4633      index[numberNonZero++]=column_[iColumn];
     4634      //index[numberNonZero++]=jColumn;
     4635    }
     4636  }
     4637  for (int iBlock=0;iBlock<numberBlocks_;iBlock++) {
     4638    // B) Can sort so just do nonbasic (and nonfixed)
     4639    // C) Can do two at a time (if so put odd one into start_)
     4640    // D) can use switch
     4641    blockStruct * block = block_ + iBlock;
     4642    //int numberPrice = block->numberInBlock_;
     4643    int numberPrice = block->numberPrice_;
     4644    int nel = block->numberElements_;
     4645    int * row = row_ + block->startElements_;
     4646    double * element = element_ + block->startElements_;
     4647    int * column = column_ + block->startIndices_;
     4648#if 0
     4649    // two at a time
     4650    if ((numberPrice&1)!=0) {
     4651      double value = 0.0;
     4652      int nel2=nel;
     4653      for (; nel2;nel2--) {
     4654        int iRow = *row++;
     4655        value += pi[iRow]*(*element++);
     4656      }
     4657      if (fabs(value)>zeroTolerance) {
     4658        array[numberNonZero]=value;
     4659        index[numberNonZero++]=*column;
     4660      }
     4661      column++;
     4662    }
     4663    numberPrice = numberPrice>>1;
     4664    switch ((nel%2)) {
     4665      // 2 k +0
     4666    case 0:
     4667      for (;numberPrice;numberPrice--) {
     4668        double value0 = 0.0;
     4669        double value1 = 0.0;
     4670        int nel2=nel;
     4671        for (; nel2;nel2--) {
     4672          int iRow0 = *row;
     4673          int iRow1 = *(row+nel);
     4674          row++;
     4675          double element0 = *element;
     4676          double element1 = *(element+nel);
     4677          element++;
     4678          value0 += pi[iRow0]*element0;
     4679          value1 += pi[iRow1]*element1;
     4680        }
     4681        row += nel;
     4682        element += nel;
     4683        if (fabs(value0)>zeroTolerance) {
     4684          array[numberNonZero]=value0;
     4685          index[numberNonZero++]=*column;
     4686        }
     4687        column++;
     4688        if (fabs(value1)>zeroTolerance) {
     4689          array[numberNonZero]=value1;
     4690          index[numberNonZero++]=*column;
     4691        }
     4692        column++;
     4693      } 
     4694      break;
     4695      // 2 k +1
     4696    case 1:
     4697      for (;numberPrice;numberPrice--) {
     4698        double value0;
     4699        double value1;
     4700        int nel2=nel-1;
     4701        {
     4702          int iRow0 = row[0];
     4703          int iRow1 = row[nel];
     4704          double pi0 = pi[iRow0];
     4705          double pi1 = pi[iRow1];
     4706          value0 = pi0*element[0];
     4707          value1 = pi1*element[nel];
     4708          row++;
     4709          element++;
     4710        }
     4711        for (; nel2;nel2--) {
     4712          int iRow0 = *row;
     4713          int iRow1 = *(row+nel);
     4714          row++;
     4715          double element0 = *element;
     4716          double element1 = *(element+nel);
     4717          element++;
     4718          value0 += pi[iRow0]*element0;
     4719          value1 += pi[iRow1]*element1;
     4720        }
     4721        row += nel;
     4722        element += nel;
     4723        if (fabs(value0)>zeroTolerance) {
     4724          array[numberNonZero]=value0;
     4725          index[numberNonZero++]=*column;
     4726        }
     4727        column++;
     4728        if (fabs(value1)>zeroTolerance) {
     4729          array[numberNonZero]=value1;
     4730          index[numberNonZero++]=*column;
     4731        }
     4732        column++;
     4733      }
     4734      break;
     4735    }
     4736#else
     4737    for (;numberPrice;numberPrice--) {
     4738      double value = 0.0;
     4739      int nel2=nel;
     4740      for (; nel2;nel2--) {
     4741        int iRow = *row++;
     4742        value += pi[iRow]*(*element++);
     4743      }
     4744      if (fabs(value)>zeroTolerance) {
     4745        array[numberNonZero]=value;
     4746        index[numberNonZero++]=*column;
     4747      }
     4748      column++;
     4749    }
     4750#endif
     4751  }
     4752  output->setNumElements(numberNonZero);
     4753}
     4754// Updates two arrays for steepest
     4755void
     4756ClpPackedMatrix3::transposeTimes2(const ClpSimplex * model,
     4757                                  const double * pi, CoinIndexedVector * output,
     4758                                  const double * piWeight,
     4759                                  double referenceIn, double devex,
     4760                                  // Array for exact devex to say what is in reference framework
     4761                                  unsigned int * reference,
     4762                                  double * weights, double scaleFactor)
     4763{
     4764  int numberNonZero=0;
     4765  int * index = output->getIndices();
     4766  double * array = output->denseVector();
     4767  double zeroTolerance = model->factorization()->zeroTolerance();
     4768  double value = 0.0;
     4769  bool killDjs = (scaleFactor==0.0);
     4770  if (!scaleFactor)
     4771    scaleFactor=1.0;
     4772  int numberOdd = block_->startIndices_;
     4773  int iColumn;
     4774  CoinBigIndex end = start_[0];
     4775  for (iColumn=0;iColumn<numberOdd;iColumn++) {
     4776    CoinBigIndex start = end;
     4777    CoinBigIndex j;
     4778    int jColumn = column_[iColumn];
     4779    end = start_[iColumn+1];
     4780    value = 0.0;
     4781    if (model->getColumnStatus(jColumn)!=ClpSimplex::basic) {
     4782      for (j=start; j<end;j++) {
     4783        int iRow = row_[j];
     4784        value -= pi[iRow]*element_[j];
     4785      }
     4786      if (fabs(value)>zeroTolerance) {
     4787        // and do other array
     4788        double modification = 0.0;
     4789        for (j=start; j<end;j++) {
     4790          int iRow = row_[j];
     4791          modification += piWeight[iRow]*element_[j];
     4792        }
     4793        double thisWeight = weights[jColumn];
     4794        double pivot = value*scaleFactor;
     4795        double pivotSquared = pivot * pivot;
     4796        thisWeight += pivotSquared * devex + pivot * modification;
     4797        if (thisWeight<DEVEX_TRY_NORM) {
     4798          if (referenceIn<0.0) {
     4799            // steepest
     4800            thisWeight = CoinMax(DEVEX_TRY_NORM,DEVEX_ADD_ONE+pivotSquared);
     4801          } else {
     4802            // exact
     4803            thisWeight = referenceIn*pivotSquared;
     4804            if (reference(jColumn))
     4805              thisWeight += 1.0;
     4806            thisWeight = CoinMax(thisWeight,DEVEX_TRY_NORM);
     4807          }
     4808        }
     4809        weights[jColumn] = thisWeight;
     4810        if (!killDjs) {
     4811          array[numberNonZero]=value;
     4812          index[numberNonZero++]=jColumn;
     4813        }
     4814      }
     4815    }
     4816  }
     4817  for (int iBlock=0;iBlock<numberBlocks_;iBlock++) {
     4818    // B) Can sort so just do nonbasic (and nonfixed)
     4819    // C) Can do two at a time (if so put odd one into start_)
     4820    // D) can use switch
     4821    blockStruct * block = block_ + iBlock;
     4822    //int numberPrice = block->numberInBlock_;
     4823    int numberPrice = block->numberPrice_;
     4824    int nel = block->numberElements_;
     4825    int * row = row_ + block->startElements_;
     4826    double * element = element_ + block->startElements_;
     4827    int * column = column_ + block->startIndices_;
     4828    for (;numberPrice;numberPrice--) {
     4829      double value = 0.0;
     4830      int nel2=nel;
     4831      for (; nel2;nel2--) {
     4832        int iRow = *row++;
     4833        value -= pi[iRow]*(*element++);
     4834      }
     4835      if (fabs(value)>zeroTolerance) {
     4836        int jColumn = *column;
     4837        // back to beginning
     4838        row -= nel;
     4839        element -= nel;
     4840        // and do other array
     4841        double modification = 0.0;
     4842        nel2=nel;
     4843        for (; nel2;nel2--) {
     4844          int iRow = *row++;
     4845          modification += piWeight[iRow]*(*element++);
     4846        }
     4847        double thisWeight = weights[jColumn];
     4848        double pivot = value*scaleFactor;
     4849        double pivotSquared = pivot * pivot;
     4850        thisWeight += pivotSquared * devex + pivot * modification;
     4851        if (thisWeight<DEVEX_TRY_NORM) {
     4852          if (referenceIn<0.0) {
     4853            // steepest
     4854            thisWeight = CoinMax(DEVEX_TRY_NORM,DEVEX_ADD_ONE+pivotSquared);
     4855          } else {
     4856            // exact
     4857            thisWeight = referenceIn*pivotSquared;
     4858            if (reference(jColumn))
     4859              thisWeight += 1.0;
     4860            thisWeight = CoinMax(thisWeight,DEVEX_TRY_NORM);
     4861          }
     4862        }
     4863        weights[jColumn] = thisWeight;
     4864        if (!killDjs) {
     4865          array[numberNonZero]=value;
     4866          index[numberNonZero++]=jColumn;
     4867        }
     4868      }
     4869      column++;
     4870    }
     4871  }
     4872  output->setNumElements(numberNonZero);
     4873  output->setPackedMode(true);
     4874}
    37194875#ifdef CLP_ALL_ONE_FILE
    37204876#undef reference
  • branches/devel/Clp/src/ClpPackedMatrix.hpp

    r800 r854  
    44#define ClpPackedMatrix_H
    55
    6 
    76#include "CoinPragma.hpp"
    87
     
    1615
    1716class ClpPackedMatrix2;
     17class ClpPackedMatrix3;
    1818class ClpPackedMatrix : public ClpMatrixBase {
    1919 
     
    260260  inline void setMatrixNull()
    261261  { matrix_=NULL;};
     262  /// Say we want special column copy
     263  inline void makeSpecialColumnCopy()
     264  { flags_ |= 8;};
     265  /// Say we don't want special column copy
     266  void releaseSpecialColumnCopy();
     267  /// Are there zeros?
     268  inline bool zeros() const
     269  { return ((flags_&1)!=0);};
     270  /// Do we want special column copy
     271  inline bool wantsSpecialColumnCopy() const
     272  { return ((flags_&8)!=0);};
    262273   //@}
    263274
     
    299310  /// make special row copy
    300311  void specialRowCopy(ClpSimplex * model,const ClpMatrixBase * rowCopy);
    301    //@}
     312  /// make special column copy
     313  void specialColumnCopy(ClpSimplex * model);
     314  /// Correct sequence in and out to give true value
     315  virtual void correctSequence(const ClpSimplex * model,int & sequenceIn, int & sequenceOut) ;
     316   //@}
     317private:
     318  /// Meat of transposeTimes by column when not scaled
     319  void gutsOfTransposeTimesUnscaled(const double * pi,CoinIndexedVector * output, const double tolerance) const;
     320  /// Meat of transposeTimes by column when scaled
     321  void gutsOfTransposeTimesScaled(const double * pi,const double * columnScale, CoinIndexedVector * output, const double tolerance) const;
     322  /// Meat of transposeTimes by row n > 2 if packed
     323  void gutsOfTransposeTimesByRowGE3(const CoinIndexedVector * piVector, CoinIndexedVector * output,
     324                                   CoinIndexedVector * spareVector, const double tolerance, const double scalar) const;
     325  /// Meat of transposeTimes by row n == 2 if packed
     326  void gutsOfTransposeTimesByRowEQ2(const CoinIndexedVector * piVector, CoinIndexedVector * output,
     327                                   CoinIndexedVector * spareVector, const double tolerance, const double scalar) const;
     328  /// Meat of transposeTimes by row n == 1 if packed
     329  void gutsOfTransposeTimesByRowEQ1(const CoinIndexedVector * piVector, CoinIndexedVector * output,
     330                                   const double tolerance, const double scalar) const;
     331  /// Gets rid of special copies
     332  void clearCopies();
    302333   
    303334   
    304335protected:
     336  /// Check validity
     337  void checkFlags() const;
    305338   /**@name Data members
    306339      The data members are protected to allow access for derived classes. */
     
    310343  /// number of active columns (normally same as number of columns)
    311344  int numberActiveColumns_;
    312   /// Zero element flag - set true if any zero elements
    313   bool zeroElements_;
    314   /// Gaps flag - set true if column start and length don't say contiguous
    315   bool hasGaps_;
     345  /** Flags -
     346      1 - has zero elements
     347      2 - has gaps
     348      4 - has special row copy
     349      8 - has special column copy
     350  */
     351  int flags_;
    316352  /// Special row copy
    317353  ClpPackedMatrix2 * rowCopy_;
     354  /// Special column copy
     355  ClpPackedMatrix3 * columnCopy_;
    318356   //@}
    319357};
     
    404442   //@}
    405443};
     444typedef struct {
     445  CoinBigIndex startElements_; // point to data
     446  int startIndices_; // point to column_
     447  int numberInBlock_;
     448  int numberPrice_; // at beginning
     449  int numberElements_; // number elements per column
     450} blockStruct;
     451class ClpPackedMatrix3 {
     452 
     453public:
     454  /**@name Useful methods */
     455  //@{
     456    /** Return <code>x * -1 * A in <code>z</code>.
     457        Note - x packed and z will be packed mode
     458        Squashes small elements and knows about ClpSimplex */
     459  void transposeTimes(const ClpSimplex * model,
     460                      const double * pi,
     461                      CoinIndexedVector * output) const;
     462  /// Updates two arrays for steepest
     463  void transposeTimes2(const ClpSimplex * model,
     464                       const double * pi, CoinIndexedVector * dj1,
     465                       const double * piWeight,
     466                       double referenceIn, double devex,
     467                       // Array for exact devex to say what is in reference framework
     468                       unsigned int * reference,
     469                       double * weights, double scaleFactor);
     470  //@}
     471
     472
     473  /**@name Constructors, destructor */
     474   //@{
     475   /** Default constructor. */
     476   ClpPackedMatrix3();
     477   /** Constructor from copy. */
     478   ClpPackedMatrix3(ClpSimplex * model,const CoinPackedMatrix * columnCopy);
     479   /** Destructor */
     480   virtual ~ClpPackedMatrix3();
     481   //@}
     482
     483   /**@name Copy method */
     484   //@{
     485   /** The copy constructor. */
     486   ClpPackedMatrix3(const ClpPackedMatrix3&);
     487   ClpPackedMatrix3& operator=(const ClpPackedMatrix3&);
     488   //@}
     489   /**@name Sort methods */
     490   //@{
     491   /** Sort blocks */
     492  void sortBlocks(const ClpSimplex * model);
     493  /// Swap one variable
     494  void swapOne(const ClpSimplex * model,const ClpPackedMatrix * matrix,
     495                int iColumn);
     496   //@}
     497   
     498   
     499protected:
     500   /**@name Data members
     501      The data members are protected to allow access for derived classes. */
     502   //@{
     503  /// Number of blocks
     504  int numberBlocks_;
     505  /// Number of columns
     506  int numberColumns_;
     507  /// Column indices and reverse lookup (within block)
     508  int * column_;
     509  /// Starts for odd/long vectors
     510  CoinBigIndex * start_;
     511  /// Rows
     512  int * row_;
     513  /// Elements
     514  double * element_;
     515  /// Blocks (ordinary start at 0 and go to first block)
     516  blockStruct * block_;
     517   //@}
     518};
    406519
    407520#endif
  • branches/devel/Clp/src/ClpPrimalColumnSteepest.cpp

    r800 r854  
    983983  if (sequenceOut>=0)
    984984    outgoingWeight=weights_[sequenceOut];
    985    
    986   // put row of tableau in rowArray and columnArray (packed)
    987   // get subset which have nonzero tableau elements
    988   transposeTimes2(updates,spareColumn1,alternateWeights_,spareColumn2,spareRow2,
    989                   -scaleFactor);
     985  // update row weights here so we can scale alternateWeights_
    990986  // update weights
    991987  double * weight;
     
    10081004    int iSequence = index[j];
    10091005    double value2 = updateBy[j];
    1010     updateBy[j]=0.0;
    10111006    ClpSimplex::Status status = model_->getStatus(iSequence+addSequence);
    10121007    double value;
     
    11221117    }
    11231118  }
     1119  // put row of tableau in rowArray and columnArray (packed)
     1120  // get subset which have nonzero tableau elements
     1121  transposeTimes2(updates,spareColumn1,alternateWeights_,spareColumn2,spareRow2,
     1122                  -scaleFactor);
     1123  // zero updateBy
     1124  CoinZeroN(updateBy,number);
    11241125  alternateWeights_->clear();
    11251126  // columns
     
    15371538    bool needSubset =  (mode_<4||numberSwitched_>1||mode_>=10);
    15381539
    1539     if (needSubset) {
    1540       // now update weight update array
    1541       model_->factorization()->updateColumnTranspose(spareRow2,
    1542                                                      alternateWeights_);
    1543       transposeTimes2(updates,spareColumn1,alternateWeights_,spareColumn2,spareRow2,0.0);
    1544     } else {
    1545       // put row of tableau in rowArray and columnArray
    1546       model_->clpMatrix()->transposeTimes(model_,-1.0,
    1547                                           updates,spareColumn2,spareColumn1);
    1548     }
    15491540    double * weight;
    15501541    double * other = alternateWeights_->denseVector();
     
    15551546    updateBy = updates->denseVector();
    15561547    weight = weights_+numberColumns;
    1557    
    15581548    if (needSubset) {
    1559       // Exact
     1549      // now update weight update array
     1550      model_->factorization()->updateColumnTranspose(spareRow2,
     1551                                                     alternateWeights_);
     1552      // do alternateWeights_ here so can scale
    15601553      for (j=0;j<number;j++) {
    15611554        int iSequence = index[j];
     
    15631556        // row has -1
    15641557        double pivot = - updateBy[j];
    1565         updateBy[j]=0.0;
    15661558        double modification = other[iSequence];
    15671559        double pivotSquared = pivot * pivot;
     
    15821574        weight[iSequence] = thisWeight;
    15831575      }
     1576      transposeTimes2(updates,spareColumn1,alternateWeights_,spareColumn2,spareRow2,0.0);
     1577    } else {
     1578      // put row of tableau in rowArray and columnArray
     1579      model_->clpMatrix()->transposeTimes(model_,-1.0,
     1580                                          updates,spareColumn2,spareColumn1);
     1581    }
     1582   
     1583    if (needSubset) {
     1584      CoinZeroN(updateBy,number);
    15841585    } else if (mode_==4) {
    15851586      // Devex
     
    27952796        alternateWeights_->capacity()==numberRows+
    27962797          model_->factorization()->maximumPivots()) {
     2798        alternateWeights_->clear();
    27972799        if (pivotSequence_>=0) {
    27982800          // save pivot order
     
    28012803          // change from pivot row number to sequence number
    28022804          pivotSequence_=pivotVariable[pivotSequence_];
    2803         } else {
    2804           alternateWeights_->clear();
    28052805        }
    28062806        state_=1;
  • branches/devel/Clp/src/ClpSimplex.cpp

    r809 r854  
    101101  savedSolution_(NULL),
    102102  numberTimesOptimal_(0),
     103  disasterArea_(NULL),
    103104  changeMade_(1),
    104105  algorithm_(0),
     
    212213  savedSolution_(NULL),
    213214  numberTimesOptimal_(0),
     215  disasterArea_(NULL),
    214216  changeMade_(1),
    215217  algorithm_(0),
     
    430432      if (!numberBasic) {
    431433        //printf("no errors on basic - going to all slack - numberOut %d\n",numberOut);
    432         allSlackBasis();
     434        allSlackBasis(true);
    433435      }
    434436      CoinSort_2(save, save + numberOut, sort,
     
    15151517  int in = sequenceIn_;
    15161518  int out = sequenceOut_;
    1517   matrix_->correctSequence(in,out);
     1519  matrix_->correctSequence(this,in,out);
    15181520  int cycle=progress_->cycle(in,out,
    15191521                            directionIn_,directionOut_);
     
    16531655  savedSolution_(NULL),
    16541656  numberTimesOptimal_(0),
     1657  disasterArea_(NULL),
    16551658  changeMade_(1),
    16561659  algorithm_(0),
     
    17581761  savedSolution_(NULL),
    17591762  numberTimesOptimal_(0),
     1763  disasterArea_(NULL),
    17601764  changeMade_(1),
    17611765  algorithm_(0),
     
    19321936  primalColumnPivot_ = rhs.primalColumnPivot_->clone(true);
    19331937  numberTimesOptimal_ = rhs.numberTimesOptimal_;
     1938  disasterArea_ = NULL;
    19341939  changeMade_ = rhs.changeMade_;
    19351940  algorithm_ = rhs.algorithm_;
     
    20552060  relaxedTolerance = relaxedTolerance +  error;
    20562061  sumOfRelaxedPrimalInfeasibilities_ = 0.0;
    2057 
    20582062  for (iRow=0;iRow<numberRows_;iRow++) {
    20592063    //assert (fabs(solution[iRow])<1.0e15||getRowStatus(iRow) == basic);
     
    25392543  if (auxiliaryModel_) {
    25402544    if (auxiliaryModel_->numberRows_==numberRows_&&
    2541         auxiliaryModel_->numberColumns_==numberColumns_)
     2545        auxiliaryModel_->numberColumns_==numberColumns_&&
     2546        (whatsChanged_&511)==511)
    25422547      oldMatrix=true;
    25432548    else
     
    27532758      if (clpMatrix&&numberThreads_)
    27542759        clpMatrix->specialRowCopy(this,rowCopy_);
     2760      if (clpMatrix)
     2761        clpMatrix->specialColumnCopy(this);
    27552762    }
    27562763  }
     
    57075714        columnActivity_[i]=0.0;
    57085715        setColumnStatus(i,atUpperBound);
     5716      }
     5717    }
     5718    if (solution_) {
     5719      // do that as well
     5720      if (!columnScale_) {
     5721        for (i=0;i<numberColumns_;i++) {
     5722          solution_[i] = columnActivity_[i];
     5723        }
     5724      } else {
     5725        for (i=0;i<numberColumns_;i++) {
     5726          solution_[i] = columnActivity_[i]*(rhsScale_/columnScale_[i]);
     5727        }
    57095728      }
    57105729    }
     
    66146633      rowArray_[0]->setNumElements(0);
    66156634      // check incoming
    6616       assert (fabs(dj_[sequenceIn_])<1.0e-6);
     6635      assert (fabs(dj_[sequenceIn_])<1.0e-6||CoinAbs(solveType_)==2);
    66176636    }
    66186637   
     
    68716890    // number of times we have declared optimality
    68726891    numberTimesOptimal_=0;
     6892    if (disasterArea_)
     6893      disasterArea_->intoSimplex();
    68736894
    68746895    return 0;
     
    70427063  otherModel.numberPrimalInfeasibilities_ = numberPrimalInfeasibilities_;
    70437064  otherModel.numberTimesOptimal_ = numberTimesOptimal_;
     7065  otherModel.disasterArea_ = NULL;
    70447066  otherModel.sumOfRelaxedDualInfeasibilities_ = sumOfRelaxedDualInfeasibilities_;
    70457067  otherModel.sumOfRelaxedPrimalInfeasibilities_ = sumOfRelaxedPrimalInfeasibilities_;
     
    71087130    iterationNumber_[i]=-1;
    71097131  }
     7132#ifdef CLP_PROGRESS_WEIGHT
     7133  for (i=0;i<CLP_PROGRESS_WEIGHT;i++) {
     7134    objectiveWeight_[i] = COIN_DBL_MAX;
     7135    infeasibilityWeight_[i] = -1.0; // set to an impossible value
     7136    realInfeasibilityWeight_[i] = COIN_DBL_MAX;
     7137    numberInfeasibilitiesWeight_[i]=-1;
     7138    iterationNumberWeight_[i]=-1;
     7139  }
     7140  drop_ =0.0;
     7141  best_ =0.0;
     7142#endif
     7143  initialWeight_=0.0;
    71107144  for (i=0;i<CLP_CYCLE;i++) {
    71117145    //obj_[i]=COIN_DBL_MAX;
     
    71377171    iterationNumber_[i]=rhs.iterationNumber_[i];
    71387172  }
     7173#ifdef CLP_PROGRESS_WEIGHT
     7174  for (i=0;i<CLP_PROGRESS_WEIGHT;i++) {
     7175    objectiveWeight_[i] = rhs.objectiveWeight_[i];
     7176    infeasibilityWeight_[i] = rhs.infeasibilityWeight_[i];
     7177    realInfeasibilityWeight_[i] = rhs.realInfeasibilityWeight_[i];
     7178    numberInfeasibilitiesWeight_[i]=rhs.numberInfeasibilitiesWeight_[i];
     7179    iterationNumberWeight_[i]=rhs.iterationNumberWeight_[i];
     7180  }
     7181  drop_ = rhs.drop_;
     7182  best_ = rhs.best_;
     7183#endif
     7184  initialWeight_ = rhs.initialWeight_;
    71397185  for (i=0;i<CLP_CYCLE;i++) {
    71407186    //obj_[i]=rhs.obj_[i];
     
    71637209    iterationNumber_[i]=-1;
    71647210  }
     7211#ifdef CLP_PROGRESS_WEIGHT
     7212  for (i=0;i<CLP_PROGRESS_WEIGHT;i++) {
     7213    objectiveWeight_[i] = COIN_DBL_MAX;
     7214    infeasibilityWeight_[i] = -1.0; // set to an impossible value
     7215    realInfeasibilityWeight_[i] = COIN_DBL_MAX;
     7216    numberInfeasibilitiesWeight_[i]=-1;
     7217    iterationNumberWeight_[i]=-1;
     7218  }
     7219  drop_ =0.0;
     7220  best_ =0.0;
     7221#endif
     7222  initialWeight_=0.0;
    71657223  for (i=0;i<CLP_CYCLE;i++) {
    71667224    //obj_[i]=COIN_DBL_MAX;
     
    71867244      iterationNumber_[i]=rhs.iterationNumber_[i];
    71877245    }
     7246#ifdef CLP_PROGRESS_WEIGHT
     7247    for (i=0;i<CLP_PROGRESS_WEIGHT;i++) {
     7248      objectiveWeight_[i] = rhs.objectiveWeight_[i];
     7249      infeasibilityWeight_[i] = rhs.infeasibilityWeight_[i];
     7250      realInfeasibilityWeight_[i] = rhs.realInfeasibilityWeight_[i];
     7251      numberInfeasibilitiesWeight_[i]=rhs.numberInfeasibilitiesWeight_[i];
     7252      iterationNumberWeight_[i]=rhs.iterationNumberWeight_[i];
     7253    }
     7254    drop_ = rhs.drop_;
     7255    best_ = rhs.best_;
     7256#endif
     7257    initialWeight_ = rhs.initialWeight_;
    71887258    for (i=0;i<CLP_CYCLE;i++) {
    71897259      //obj_[i]=rhs.obj_[i];
     
    87758845  //}
    87768846}
     8847// Compute minimization objective value from internal solution
     8848double
     8849ClpSimplex::computeInternalObjectiveValue()
     8850{
     8851  int iSequence;
     8852  //double oldObj = objectiveValue_;
     8853  double objectiveValue = -dblParam_[ClpObjOffset];
     8854  const double * obj = objective();
     8855  double scaleR = optimizationDirection_/rhsScale_;
     8856  if (!columnScale_) {
     8857    for (iSequence=0;iSequence<numberColumns_;iSequence++) {
     8858      double value = solution_[iSequence]*scaleR;
     8859      objectiveValue += value*obj[iSequence];
     8860    }
     8861  } else {
     8862    for (iSequence=0;iSequence<numberColumns_;iSequence++) {
     8863      double value = solution_[iSequence]*scaleR*columnScale_[iSequence];
     8864      objectiveValue += value*obj[iSequence];
     8865    }
     8866  }
     8867  return objectiveValue;
     8868}
     8869void
     8870ClpSimplex::setSpecialOptions(unsigned int value)
     8871{
     8872  specialOptions_=value;
     8873}
    87778874// If user left factorization frequency then compute
    87788875void
  • branches/devel/Clp/src/ClpSimplex.hpp

    r805 r854  
    1010#ifndef ClpSimplex_H
    1111#define ClpSimplex_H
    12 
    13 #if defined(_MSC_VER)
    14 // Turn off compiler warning about long names
    15 #  pragma warning(disable:4786)
    16 #endif
    1712
    1813#include <iostream>
     
    3025class OsiClpSolverInterface;
    3126class CoinWarmStartBasis;
     27class ClpDisasterHandler;
    3228
    3329/** This solves LPs using the simplex method
     
    478474  inline void setNumberDualInfeasibilities(int value)
    479475          { numberDualInfeasibilities_=value;} ;
     476  /// Number of dual infeasibilities (without free)
     477  inline int numberDualInfeasibilitiesWithoutFree() const
     478          { return numberDualInfeasibilitiesWithoutFree_;} ;
    480479  /// Sum of primal infeasibilities
    481480  inline double sumPrimalInfeasibilities() const
     
    664663          { return largestSolutionError_;} ;
    665664public:
     665  /// Disaster handler
     666  inline void setDisasterHandler(ClpDisasterHandler * handler)
     667  { disasterArea_= handler;};
    666668  /// Large bound value (for complementarity etc)
    667669  inline double largeValue() const
     
    928930   /// Compute objective value from solution and put in objectiveValue_
    929931  void computeObjectiveValue();
     932  /// Compute minimization objective value from internal solution without perturbation
     933  double computeInternalObjectiveValue();
    930934  /** Number of extra rows.  These are ones which will be dynamically created
    931935      each iteration.  This is for GUB but may have other uses.
     
    964968             repository.  See COIN_CLP_VETTED comments.
    965969      0x01000000 is Cbc (and in branch and bound)
     970      0x02000000 is in a different branch and bound
    966971  */
    967972#define COIN_CBC_USING_CLP 0x01000000
    968973  inline unsigned int specialOptions() const
    969974  { return specialOptions_;};
    970   inline void setSpecialOptions(unsigned int value)
    971   { specialOptions_=value;};
     975  void setSpecialOptions(unsigned int value);
    972976  //@}
    973977
     
    12201224  /// Number of times code has tentatively thought optimal
    12211225  int numberTimesOptimal_;
     1226  /// Disaster handler
     1227  ClpDisasterHandler * disasterArea_;
    12221228  /// If change has been made (first attempt at stopping looping)
    12231229  int changeMade_;
     
    13681374  /**@name Data  */
    13691375#define CLP_PROGRESS 5
     1376  //#define CLP_PROGRESS_WEIGHT 10
    13701377  //@{
    13711378  /// Objective values
     
    13751382  /// Sum of real primal infeasibilities for primal
    13761383  double realInfeasibility_[CLP_PROGRESS];
     1384#ifdef CLP_PROGRESS_WEIGHT
     1385  /// Objective values for weights
     1386  double objectiveWeight_[CLP_PROGRESS_WEIGHT];
     1387  /// Sum of infeasibilities for algorithm for weights
     1388  double infeasibilityWeight_[CLP_PROGRESS_WEIGHT];
     1389  /// Sum of real primal infeasibilities for primal for weights
     1390  double realInfeasibilityWeight_[CLP_PROGRESS_WEIGHT];
     1391  /// Drop  for weights
     1392  double drop_;
     1393  /// Best? for weights
     1394  double best_;
     1395#endif
     1396  /// Initial weight for weights
     1397  double initialWeight_;
    13771398#define CLP_CYCLE 12
    13781399  /// For cycle checking
     
    13871408  /// Iteration number at which occurred
    13881409  int iterationNumber_[CLP_PROGRESS];
     1410#ifdef CLP_PROGRESS_WEIGHT
     1411  /// Number of infeasibilities for weights
     1412  int numberInfeasibilitiesWeight_[CLP_PROGRESS_WEIGHT];
     1413  /// Iteration number at which occurred for weights
     1414  int iterationNumberWeight_[CLP_PROGRESS_WEIGHT];
     1415#endif
    13891416  /// Number of times checked (so won't stop too early)
    13901417  int numberTimes_;
  • branches/devel/Clp/src/ClpSimplexDual.cpp

    r822 r854  
    384384      }
    385385    }
     386    // see if in Cbc etc
     387    bool inCbcOrOther = (specialOptions_&0x03000000)!=0;
     388#if 0
     389    bool gotoPrimal=false;
     390    if (inCbcOrOther&&numberIterations_>disasterArea_+numberRows_&&
     391        numberDualInfeasibilitiesWithoutFree_&&largestDualError_>1.0e-1) {
     392      if (!disasterArea_) {
     393        printf("trying all slack\n");
     394        // try all slack basis
     395        allSlackBasis(true);
     396        disasterArea_=2*numberRows_;
     397      } else {
     398        printf("going to primal\n");
     399        // go to primal
     400        gotoPrimal=true;
     401        allSlackBasis(true);
     402      }
     403    }
     404#endif
     405    bool disaster=false;
     406    if (disasterArea_&&inCbcOrOther&&disasterArea_->check()) {
     407      disasterArea_->saveInfo();
     408      disaster=true;
     409    }
    386410    // may factorize, checks if problem finished
    387411    statusOfProblemInDual(lastCleaned,factorType,saveDuals,data,
    388412                          ifValuesPass);
     413    if (disaster)
     414      problemStatus_=3;
    389415    // If values pass then do easy ones on first time
    390416    if (ifValuesPass&&
     
    900926#if CLP_DEBUG>2
    901927    // very expensive
    902     if (numberIterations_>0&&numberIterations_<-801) {
    903       handler_->setLogLevel(63);
     928    if (numberIterations_>3063&&numberIterations_<30700) {
     929      //handler_->setLogLevel(63);
    904930      double saveValue = objectiveValue_;
    905931      double * saveRow1 = new double[numberRows_];
     
    14181444          printf("obj backwards %g %g\n",objectiveValue_,oldobj);
    14191445#endif
     1446#if 0
     1447        {
     1448          for (int i=0;i<numberRows_+numberColumns_;i++) {
     1449            FakeBound bound = getFakeBound(i);
     1450            if (bound==ClpSimplexDual::upperFake) {
     1451              assert (upper_[i]<1.0e20);
     1452            } else if (bound==ClpSimplexDual::lowerFake) {
     1453              assert (lower_[i]>-1.0e20);
     1454            } else if (bound==ClpSimplexDual::bothFake) {
     1455              assert (upper_[i]<1.0e20);
     1456              assert (lower_[i]>-1.0e20);
     1457            }
     1458          }
     1459        }
     1460#endif
    14201461        if (whatNext==1||candidate==-2) {
    14211462          problemStatus_ =-2; // refactorize
     
    23992440                           CoinIndexedVector * spareArray,
    24002441                           double acceptablePivot,
    2401                            double & upperReturn, double &bestReturn)
     2442                           double & upperReturn, double &bestReturn,double & badFree)
    24022443{
    24032444  // do first pass to get possibles
     
    24152456  int numberRemaining=0;
    24162457  int i;
     2458  badFree=0.0;
    24172459  for (int iSection=0;iSection<2;iSection++) {
    24182460
     
    24502492        bestPossible = CoinMax(bestPossible,fabs(alpha));
    24512493        oldValue = reducedCost[iSequence];
    2452         // If freehas to be very large - should come in via dualRow
    2453         if (getStatus(iSequence+addSequence)==isFree&&fabs(alpha)<1.0e-3)
    2454           break;
     2494        // If free has to be very large - should come in via dualRow
     2495        //if (getStatus(iSequence+addSequence)==isFree&&fabs(alpha)<1.0e-3)
     2496        //break;
    24552497        if (oldValue>dualTolerance_) {
    24562498          keep = true;
     
    24582500          keep = true;
    24592501        } else {
    2460           if (fabs(alpha)>CoinMax(10.0*acceptablePivot,1.0e-5))
     2502          if (fabs(alpha)>CoinMax(10.0*acceptablePivot,1.0e-5)) {
    24612503            keep = true;
    2462           else
     2504          } else {
    24632505            keep=false;
     2506            badFree=CoinMax(badFree,fabs(alpha));
     2507          }
    24642508        }
    24652509        if (keep) {
     
    25962640  upperTheta = 1.0e31;
    25972641  double bestPossible=0.0;
     2642  double badFree=0.0;
    25982643  if (spareIntArray_[0]!=-1) {
    25992644    numberRemaining = dualColumn0(rowArray,columnArray,spareArray,
    2600                                   acceptablePivot,upperTheta,bestPossible);
     2645                                  acceptablePivot,upperTheta,bestPossible,badFree);
    26012646  } else {
    26022647    // already done
     
    31093154      alpha_=0.0;
    31103155    }
     3156    if (fabs(theta_*badFree)>10.0*dualTolerance_&&factorization_->pivots()) {
     3157      if (handler_->logLevel()>1)
     3158        printf("forcing re-factorizationon free\n");
     3159      alpha_=0.0;
     3160    }
    31113161#endif
    31123162    lowerIn_ = lower_[sequenceIn_];
     
    32713321    forceFactorization_=1; // a bit drastic but ..
    32723322    changeMade_++; // say something changed
     3323    // get correct bounds on all variables
     3324    resetFakeBounds();
    32733325  }
    32743326  int tentativeStatus = problemStatus_;
     
    33033355        matrix_->generalExpanded(this,6,dummy);
    33043356        // get correct bounds on all variables
    3305         double dummyChangeCost=0.0;
    3306         changeBounds(true,rowArray_[2],dummyChangeCost);
    3307         // throw away change
    3308         for (i=0;i<4;i++)
    3309           rowArray_[i]->clear();
     3357        resetFakeBounds();
    33103358        // need to reject something
    33113359        char x = isColumn(sequenceOut_) ? 'C' :'R';
     
    33443392            factorization_->setDenseThreshold(saveDense);
    33453393          }
     3394          // get correct bounds on all variables
     3395          resetFakeBounds();
    33463396        }
    33473397      }
     
    33733423    matrix_->generalExpanded(this,6,dummy);
    33743424    // get correct bounds on all variables
    3375     //double dummyChangeCost=0.0;
    3376     //changeBounds(true,rowArray_[2],dummyChangeCost);
    3377     // throw away change
    3378     for (i=0;i<4;i++)
    3379       rowArray_[i]->clear();
     3425    resetFakeBounds();
    33803426    // need to reject something
    33813427    char x = isColumn(sequenceOut_) ? 'C' :'R';
     
    34153461        factorization_->setDenseThreshold(saveDense);
    34163462      }
     3463      // get correct bounds on all variables
     3464      resetFakeBounds();
    34173465    }
    34183466    // get primal and dual solutions
     
    34543502        matrix_->generalExpanded(this,6,dummy);
    34553503        // get correct bounds on all variables
    3456         double dummyChangeCost=0.0;
    3457         changeBounds(true,rowArray_[2],dummyChangeCost);
    3458         // throw away change
    3459         for (int i=0;i<4;i++)
    3460           rowArray_[i]->clear();
     3504        resetFakeBounds();
    34613505        if(factorization_->pivotTolerance()<0.2)
    34623506          factorization_->pivotTolerance(0.2);
     
    34843528            factorization_->setDenseThreshold(saveDense);
    34853529          }
     3530          resetFakeBounds();
    34863531        }
    34873532        type = 2; // so will restore weights
     
    35963641          <<CoinMessageEol;
    35973642        // save solution in case unbounded
    3598         CoinMemcpyN(columnActivityWork_,numberColumns_,
    3599                           columnArray_[0]->denseVector());
    3600         CoinMemcpyN(rowActivityWork_,numberRows_,
    3601                           rowArray_[2]->denseVector());
     3643        double * saveColumnSolution = NULL;
     3644        double * saveRowSolution = NULL;
     3645        bool inCbc = (specialOptions_&0x01000000)!=0;
     3646        if (!inCbc) {
     3647          saveColumnSolution = CoinCopyOfArray(columnActivityWork_,numberColumns_);
     3648          saveRowSolution = CoinCopyOfArray(rowActivityWork_,numberRows_);
     3649        }
    36023650        numberChangedBounds=changeBounds(false,rowArray_[3],changeCost);
    36033651        if (numberChangedBounds<=0&&!numberDualInfeasibilities_) {
     
    36263674            if (numberTimesOptimal_==1) {
    36273675              dualTolerance_ = dblParam_[ClpDualTolerance];
    3628               // better to have small tolerance even if slower
    3629               factorization_->zeroTolerance(1.0e-15);
    36303676            } else {
     3677              if (numberTimesOptimal_==2) {
     3678                // better to have small tolerance even if slower
     3679                factorization_->zeroTolerance(1.0e-15);
     3680              }
    36313681              dualTolerance_ = dblParam_[ClpDualTolerance];
    36323682              dualTolerance_ *= pow(2.0,numberTimesOptimal_-1);
     
    37003750                if(getColumnStatus(iColumn)!= basic)
    37013751                  ray_[iColumn] +=
    3702                     columnActivityWork_[iColumn]-original[iColumn];
     3752                    saveColumnSolution[iColumn]-original[iColumn];
    37033753                columnActivityWork_[iColumn] = original[iColumn];
    37043754              }
    3705               CoinMemcpyN(rowArray_[2]->denseVector(),numberRows_,
     3755              CoinMemcpyN(saveRowSolution,numberRows_,
    37063756                                rowActivityWork_);
    37073757            }
     
    37113761          }
    37123762        }
    3713         CoinZeroN(columnArray_[0]->denseVector(),numberColumns_);
    3714         CoinZeroN(rowArray_[2]->denseVector(),numberRows_);
     3763        delete [] saveColumnSolution;
     3764        delete [] saveRowSolution;
    37153765      }
    37163766      if (problemStatus_==-4||problemStatus_==-5) {
     
    39684018  double limit = 0.0;
    39694019  getDblParam(ClpDualObjectiveLimit, limit);
     4020#if 0
    39704021  if(fabs(limit)<1.0e30&&objectiveValue()*optimizationDirection_>
    3971            limit&&
    3972            !numberAtFakeBound()&&!numberDualInfeasibilities_) {
    3973     //printf("lim %g obj %g %g\n",limit,objectiveValue_,objectiveValue());
    3974     if (perturbation_==101) {
    3975       // be careful
    3976       if (numberIterations_) {
    3977         computeObjectiveValue(); // value without perturbation
    3978         if(objectiveValue()*optimizationDirection_>limit) {
    3979           problemStatus_=1;
    3980           secondaryStatus_ = 1; // and say was on cutoff
    3981         }
    3982       }
    3983     } else {
    3984       problemStatus_=1;
    3985       secondaryStatus_ = 1; // and say was on cutoff
     4022     limit+100.0) {
     4023    printf("lim %g obj %g %g - wo perturb %g sum dual %g\n",
     4024           limit,objectiveValue_,objectiveValue(),computeInternalObjectiveValue(),sumDualInfeasibilities_);
     4025  }
     4026#endif
     4027  if(fabs(limit)<1.0e30&&objectiveValue()*optimizationDirection_>
     4028     limit&&!numberAtFakeBound()) {
     4029    bool looksInfeasible = !numberDualInfeasibilities_;
     4030    if (objectiveValue()*optimizationDirection_>limit+fabs(0.1*limit)+1.0e2*sumDualInfeasibilities_+1.0e4&&
     4031        sumDualInfeasibilities_<largestDualError_&&numberIterations_>0.5*numberRows_+1000)
     4032      looksInfeasible=true;
     4033    if (looksInfeasible) {
     4034      //printf("lim %g obj %g %g\n",limit,objectiveValue_,objectiveValue());
     4035      if (perturbation_==101) {
     4036        // be careful
     4037        if (numberIterations_) {
     4038          if(computeInternalObjectiveValue()>limit) {
     4039            problemStatus_=1;
     4040            secondaryStatus_ = 1; // and say was on cutoff
     4041          }
     4042        }
     4043      } else {
     4044        problemStatus_=1;
     4045        secondaryStatus_ = 1; // and say was on cutoff
     4046      }
    39864047    }
    39874048  }
     
    40054066        secondaryStatus_ = 1; // and say was on cutoff
    40064067      } else if (largestPrimalError_>1.0e5) {
    4007         allSlackBasis();
     4068        {
     4069          int iBigB=-1;
     4070          double bigB=0.0;
     4071          int iBigN=-1;
     4072          double bigN=0.0;
     4073          for (int i=0;i<numberRows_+numberColumns_;i++) {
     4074            double value = fabs(solution_[i]);
     4075            if (getStatus(i)==basic) {
     4076              if (value>bigB) {
     4077                bigB= value;
     4078                iBigB=i;
     4079              }
     4080            } else {
     4081              if (value>bigN) {
     4082                bigN= value;
     4083                iBigN=i;
     4084              }
     4085            }
     4086          }
     4087          if (bigB>1.0e8||bigN>1.0e8) {
     4088            if (handler_->logLevel()>0)
     4089              printf("it %d - basic %d %g, nonbasic %d %g\n",
     4090                     numberIterations_,iBigB,bigB,iBigN,bigN);
     4091          }
     4092        }
     4093        if (handler_->logLevel()>0)
     4094          printf("bad dual - going to primal %d %g\n",looksBad,largestPrimalError_);
     4095        allSlackBasis(true);
    40084096        problemStatus_=10;
    4009         //if (handler_->logLevel()>0)
    4010           printf("bad dual - going to primal %d %g\n",looksBad,largestPrimalError_);
    40114097      }
    40124098    }
     
    40354121  }
    40364122#endif
     4123  // Allow matrices to be sorted etc
     4124  int fake=-999; // signal sort
     4125  matrix_->correctSequence(this,fake,fake);
    40374126}
    40384127/* While updateDualsInDual sees what effect is of flip
     
    41034192    lower_[iSequence]=auxiliaryModel_->lowerRegion()[iSequence+numberRows_+numberColumns_];
    41044193    upper_[iSequence]=auxiliaryModel_->upperRegion()[iSequence+numberRows_+numberColumns_];
     4194  setFakeBound(iSequence,noFake);
    41054195    return;
    41064196  }
     
    56195709  }
    56205710}
     5711void
     5712ClpSimplexDual::resetFakeBounds()
     5713{
     5714  double dummyChangeCost=0.0;
     5715  changeBounds(true,rowArray_[2],dummyChangeCost);
     5716  // throw away change
     5717  for (int i=0;i<4;i++)
     5718    rowArray_[i]->clear();
     5719}
  • branches/devel/Clp/src/ClpSimplexDual.hpp

    r754 r854  
    201201                  CoinIndexedVector * spareArray,
    202202                  double acceptablePivot,
    203                   double & upperReturn, double &bestReturn);
     203                  double & upperReturn, double &bestReturn,double & badFree);
    204204  /**
    205205      Row array has row part of pivot row
     
    293293                  ClpDataSave & saveData);
    294294  //int dual2(int ifValuesPass,int startFinishOptions=0);
     295  void resetFakeBounds();
    295296 
    296297  //@}
  • branches/devel/Clp/src/ClpSimplexNonlinear.cpp

    r800 r854  
    589589  //if (goToDual)
    590590  //problemStatus_=10; // try dual
     591  // Allow matrices to be sorted etc
     592  int fake=-999; // signal sort
     593  matrix_->correctSequence(this,fake,fake);
    591594}
    592595/*
  • branches/devel/Clp/src/ClpSimplexOther.cpp

    r754 r854  
    20022002      numberDualInfeasibilities_=1;
    20032003  }
     2004  // Allow matrices to be sorted etc
     2005  int fake=-999; // signal sort
     2006  matrix_->correctSequence(this,fake,fake);
    20042007}
    20052008/* This has the flow between re-factorizations
  • branches/devel/Clp/src/ClpSimplexPrimal.cpp

    r844 r854  
    433433  }
    434434  // if infeasible get real values
     435  //printf("XXXXY final cost %g\n",infeasibilityCost_);
     436  progress_->initialWeight_=0.0;
    435437  if (problemStatus_==1) {
    436438    infeasibilityCost_=0.0;
     
    878880    nonLinearCost_->checkInfeasibilities(primalTolerance_);
    879881  }
     882  if (nonLinearCost_->numberInfeasibilities()>0&&!progress->initialWeight_&&!ifValuesPass&&infeasibilityCost_==1.0e10) {
     883    // first time infeasible - start up weight computation
     884    double * oldDj = dj_;
     885    double * oldCost = cost_;
     886    int numberRows2 = numberRows_+numberExtraRows_;
     887    int numberTotal = numberRows2+numberColumns_;
     888    dj_ = new double[numberTotal];
     889    cost_ = new double[numberTotal];
     890    reducedCostWork_ = dj_;
     891    rowReducedCost_ = dj_+numberColumns_;
     892    objectiveWork_ = cost_;
     893    rowObjectiveWork_ = cost_+numberColumns_;
     894    double direction = optimizationDirection_*objectiveScale_;
     895    const double * obj = objective();
     896    memset(rowObjectiveWork_,0,numberRows_*sizeof(double));
     897    int iSequence;
     898    if (columnScale_)
     899      for (iSequence=0;iSequence<numberColumns_;iSequence++)
     900        cost_[iSequence] = obj[iSequence]*direction*columnScale_[iSequence];
     901    else
     902      for (iSequence=0;iSequence<numberColumns_;iSequence++)
     903        cost_[iSequence] = obj[iSequence]*direction;
     904    computeDuals(NULL);
     905    int numberSame=0;
     906    int numberDifferent=0;
     907    int numberZero=0;
     908    int numberFreeSame=0;
     909    int numberFreeDifferent=0;
     910    int numberFreeZero=0;
     911    int n=0;
     912    for (iSequence=0;iSequence<numberTotal;iSequence++) {
     913      if (getStatus(iSequence) != basic&&!flagged(iSequence)) {
     914        // not basic
     915        double distanceUp = upper_[iSequence]-solution_[iSequence];
     916        double distanceDown = solution_[iSequence]-lower_[iSequence];
     917        double feasibleDj = dj_[iSequence];
     918        double infeasibleDj = oldDj[iSequence]-feasibleDj;
     919        double value = feasibleDj*infeasibleDj;
     920        if (distanceUp>primalTolerance_) {
     921          // Check if "free"
     922          if (distanceDown>primalTolerance_) {
     923            // free
     924            if (value>dualTolerance_) {
     925              numberFreeSame++;
     926            } else if(value<-dualTolerance_) {
     927              numberFreeDifferent++;
     928              dj_[n++] = feasibleDj/infeasibleDj;
     929            } else {
     930              numberFreeZero++;
     931            }
     932          } else {
     933            // should not be negative
     934            if (value>dualTolerance_) {
     935              numberSame++;
     936            } else if(value<-dualTolerance_) {
     937              numberDifferent++;
     938              dj_[n++] = feasibleDj/infeasibleDj;
     939            } else {
     940              numberZero++;
     941            }
     942          }
     943        } else if (distanceDown>primalTolerance_) {
     944          // should not be positive
     945          if (value>dualTolerance_) {
     946              numberSame++;
     947            } else if(value<-dualTolerance_) {
     948              numberDifferent++;
     949              dj_[n++] = feasibleDj/infeasibleDj;
     950            } else {
     951              numberZero++;
     952            }
     953        }
     954      }
     955      progress->initialWeight_=-1.0;
     956    }
     957    //printf("XXXX %d same, %d different, %d zero, -- free %d %d %d\n",
     958    //   numberSame,numberDifferent,numberZero,
     959    //   numberFreeSame,numberFreeDifferent,numberFreeZero);
     960    // we want most to be same
     961    if (n) {
     962      double most = 0.95;
     963      std::sort(dj_,dj_+n);
     964      int which = (int) ((1.0-most)*((double) n));
     965      double take = -dj_[which]*infeasibilityCost_;
     966      //printf("XXXXZ inf cost %g take %g (range %g %g)\n",infeasibilityCost_,take,-dj_[0]*infeasibilityCost_,-dj_[n-1]*infeasibilityCost_);
     967      take = -dj_[0]*infeasibilityCost_;
     968      infeasibilityCost_ = CoinMin(1000.0*take,1.0000001e10);;
     969      //printf("XXXX increasing weight to %g\n",infeasibilityCost_);
     970    }
     971    delete [] dj_;
     972    delete [] cost_;
     973    dj_= oldDj;
     974    cost_ = oldCost;
     975    reducedCostWork_ = dj_;
     976    rowReducedCost_ = dj_+numberColumns_;
     977    objectiveWork_ = cost_;
     978    rowObjectiveWork_ = cost_+numberColumns_;
     979    if (n)
     980      gutsOfSolution(NULL,NULL,ifValuesPass!=0);
     981  }
    880982  double trueInfeasibility =nonLinearCost_->sumInfeasibilities();
    881   if (!nonLinearCost_->numberInfeasibilities()&&infeasibilityCost_==1.0e10&&!ifValuesPass) {
     983  if (!nonLinearCost_->numberInfeasibilities()&&infeasibilityCost_==1.0e10&&!ifValuesPass&&true) {
    882984    // relax if default
    883     infeasibilityCost_ = CoinMin(CoinMax(100.0*sumDualInfeasibilities_,1.0e4),1.0e7);
     985    infeasibilityCost_ = CoinMin(CoinMax(100.0*sumDualInfeasibilities_,1.0e6),1.00000001e10);
    884986    // reset looping criterion
    885987    *progress = ClpSimplexProgress();
     
    9151017      numberDualInfeasibilities_ = -776;
    9161018    }
    917   }
     1019    // But if real primal infeasibilities nonzero carry on
     1020    if (nonLinearCost_->numberInfeasibilities()) {
     1021      // most likely to happen if infeasible
     1022      double relaxedToleranceP=primalTolerance_;
     1023      // we can't really trust infeasibilities if there is primal error
     1024      double error = CoinMin(1.0e-2,largestPrimalError_);
     1025      // allow tolerance at least slightly bigger than standard
     1026      relaxedToleranceP = relaxedToleranceP +  error;
     1027      int ninfeas = nonLinearCost_->numberInfeasibilities();
     1028      double sum = nonLinearCost_->sumInfeasibilities();
     1029      double average = sum/ ((double) ninfeas);
     1030#ifdef COIN_DEVELOP
     1031      if (handler_->logLevel()>0)
     1032        printf("nonLinearCost says infeasible %d summing to %g\n",
     1033               ninfeas,sum);
     1034#endif
     1035      if (average>relaxedToleranceP) {
     1036        sumOfRelaxedPrimalInfeasibilities_ = sum;
     1037        numberPrimalInfeasibilities_ = ninfeas;
     1038        sumPrimalInfeasibilities_ = sum;
     1039      }
     1040    }
     1041  }
    9181042  // had ||(type==3&&problemStatus_!=-5) -- ??? why ????
    9191043  if ((dualFeasible()||problemStatus_==-4)&&!ifValuesPass) {
     
    12011325  }
    12021326#endif
     1327  // Allow matrices to be sorted etc
     1328  int fake=-999; // signal sort
     1329  matrix_->correctSequence(this,fake,fake);
    12031330}
    12041331/*
  • branches/devel/Clp/src/ClpSolve.cpp

    r809 r854  
    33
    44// This file has higher level solve functions
    5 
    65
    76#include "ClpConfig.h"
     
    965964      if (model2->numberRows()>100)
    966965        model2->setSpecialOptions(saveOptions|64); // go as far as possible
    967       model2->dual(0);
    968       model2->setSpecialOptions(saveOptions);
     966      //int numberRows = model2->numberRows();
     967      //int numberColumns = model2->numberColumns();
     968      if (dynamic_cast< ClpPackedMatrix*>(matrix_)) {
     969        // See if original wanted vector
     970        ClpPackedMatrix * clpMatrixO = dynamic_cast< ClpPackedMatrix*>(matrix_);
     971        ClpMatrixBase * matrix = model2->clpMatrix();
     972        if (dynamic_cast< ClpPackedMatrix*>(matrix)&&clpMatrixO->wantsSpecialColumnCopy()) {
     973          ClpPackedMatrix * clpMatrix = dynamic_cast< ClpPackedMatrix*>(matrix);
     974          clpMatrix->makeSpecialColumnCopy();
     975          //model2->setSpecialOptions(model2->specialOptions()|256); // to say no row copy for comparisons
     976          model2->dual(0);
     977          clpMatrix->releaseSpecialColumnCopy();
     978        } else {
     979          model2->dual(0);
     980        }
     981      } else {
     982        model2->dual(0);
     983      }
    969984    } else if (!numberNotE&&0) {
    970985      // E so we can do in another way
     
    12491264    }
    12501265#endif
    1251     model2->primal(primalStartup);
     1266    if (dynamic_cast< ClpPackedMatrix*>(matrix_)) {
     1267      // See if original wanted vector
     1268      ClpPackedMatrix * clpMatrixO = dynamic_cast< ClpPackedMatrix*>(matrix_);
     1269      ClpMatrixBase * matrix = model2->clpMatrix();
     1270      if (dynamic_cast< ClpPackedMatrix*>(matrix)&&clpMatrixO->wantsSpecialColumnCopy()) {
     1271        ClpPackedMatrix * clpMatrix = dynamic_cast< ClpPackedMatrix*>(matrix);
     1272        clpMatrix->makeSpecialColumnCopy();
     1273        //model2->setSpecialOptions(model2->specialOptions()|256); // to say no row copy for comparisons
     1274        model2->primal(primalStartup);
     1275        clpMatrix->releaseSpecialColumnCopy();
     1276      } else {
     1277        model2->primal(primalStartup);
     1278      }
     1279    } else {
     1280      model2->primal(primalStartup);
     1281    }
    12521282    time2 = CoinCpuTime();
    12531283    timeCore = time2-timeX;
     
    15451575      if (interrupt)
    15461576        currentModel = &small;
    1547       small.primal(1);
     1577      if (dynamic_cast< ClpPackedMatrix*>(matrix_)) {
     1578        // See if original wanted vector
     1579        ClpPackedMatrix * clpMatrixO = dynamic_cast< ClpPackedMatrix*>(matrix_);
     1580        ClpMatrixBase * matrix = small.clpMatrix();
     1581        if (dynamic_cast< ClpPackedMatrix*>(matrix)&&clpMatrixO->wantsSpecialColumnCopy()) {
     1582          ClpPackedMatrix * clpMatrix = dynamic_cast< ClpPackedMatrix*>(matrix);
     1583          clpMatrix->makeSpecialColumnCopy();
     1584          small.primal(1);
     1585          clpMatrix->releaseSpecialColumnCopy();
     1586        } else {
     1587          small.primal(1);
     1588        }
     1589      } else {
     1590        small.primal(1);
     1591      }
    15481592      totalIterations += small.numberIterations();
    15491593      // move solution back
     
    21482192      primal(1);
    21492193      setPerturbation(savePerturbation);
    2150       numberIterations += this->numberIterations();
     2194      numberIterations += numberIterations_;
     2195      numberIterations_ = numberIterations;
    21512196      finalStatus=status();
    21522197      time2 = CoinCpuTime();
  • branches/devel/Clp/src/Idiot.cpp

    r754 r854  
    10131013    /*printf("%d in basis\n",ninbas);*/
    10141014  }
     1015  bool wantVector=false;
     1016  if (dynamic_cast< ClpPackedMatrix*>(model_->clpMatrix())) {
     1017    // See if original wanted vector
     1018    ClpPackedMatrix * clpMatrixO = dynamic_cast< ClpPackedMatrix*>(model_->clpMatrix());
     1019    wantVector = clpMatrixO->wantsSpecialColumnCopy();
     1020  }
    10151021  if (addAll<3) {
    10161022    ClpPresolve pinfo;
     
    10201026    }
    10211027    if (model_) {
    1022       model_->primal(1);
     1028      if (!wantVector) {
     1029        model_->primal(1);
     1030      } else {
     1031        ClpMatrixBase * matrix = model_->clpMatrix();
     1032        ClpPackedMatrix * clpMatrix = dynamic_cast< ClpPackedMatrix*>(matrix);
     1033        assert (clpMatrix);
     1034        clpMatrix->makeSpecialColumnCopy();
     1035        model_->primal(1);
     1036        clpMatrix->releaseSpecialColumnCopy();
     1037      }
    10231038      if (presolve) {
    10241039        pinfo.postsolve(true);
     
    10671082        presolve=0;
    10681083      }
    1069       model_->primal(1);
     1084      if (!wantVector) {
     1085        model_->primal(1);
     1086      } else {
     1087        ClpMatrixBase * matrix = model_->clpMatrix();
     1088        ClpPackedMatrix * clpMatrix = dynamic_cast< ClpPackedMatrix*>(matrix);
     1089        assert (clpMatrix);
     1090        clpMatrix->makeSpecialColumnCopy();
     1091        model_->primal(1);
     1092        clpMatrix->releaseSpecialColumnCopy();
     1093      }
    10701094      if (presolve) {
    10711095        pinfo.postsolve(true);
     
    10921116        presolve=0;
    10931117      }
    1094       model_->primal(1);
     1118      if (!wantVector) {
     1119        model_->primal(1);
     1120      } else {
     1121        ClpMatrixBase * matrix = model_->clpMatrix();
     1122        ClpPackedMatrix * clpMatrix = dynamic_cast< ClpPackedMatrix*>(matrix);
     1123        assert (clpMatrix);
     1124        clpMatrix->makeSpecialColumnCopy();
     1125        model_->primal(1);
     1126        clpMatrix->releaseSpecialColumnCopy();
     1127      }
    10951128      if (presolve) {
    10961129        pinfo.postsolve(true);
  • branches/devel/Clp/src/unitTest.cpp

    r799 r854  
    241241//----------------------------------------------------------------
    242242int mainTest (int argc, const char *argv[],int algorithm,
    243               ClpSimplex empty, bool doPresolve, int switchOffValue)
     243              ClpSimplex empty, bool doPresolve, int switchOffValue,bool doVector)
    244244{
    245245  int i;
     
    503503            if (solution.maximumSeconds()<0.0)
    504504              solution.setMaximumSeconds(120.0);
     505            if (doVector) {
     506              ClpMatrixBase * matrix = solution.clpMatrix();
     507              if (dynamic_cast< ClpPackedMatrix*>(matrix)) {
     508                ClpPackedMatrix * clpMatrix = dynamic_cast< ClpPackedMatrix*>(matrix);
     509                clpMatrix->makeSpecialColumnCopy();
     510              }
     511            }
    505512            solution.initialSolve(solveOptions);
    506513            double time2 = CoinCpuTime()-time1;
     
    591598          solveOptions.setPresolveType(ClpSolve::presolveOff);
    592599      }
     600      if (doVector) {
     601        ClpMatrixBase * matrix = solution.clpMatrix();
     602        if (dynamic_cast< ClpPackedMatrix*>(matrix)) {
     603          ClpPackedMatrix * clpMatrix = dynamic_cast< ClpPackedMatrix*>(matrix);
     604          clpMatrix->makeSpecialColumnCopy();
     605        }
     606      }
    593607      solution.initialSolve(solveOptions);
    594608      double time2 = CoinCpuTime()-time1;
     
    784798    model.setFactorizationFrequency(10);
    785799    model.primal();
    786 
     800    model.primal(0,3);
     801    model.setObjCoeff(3,-2.9473684210526314);
     802    model.primal(0,3);
    787803    // Write saved solutions
    788804    int nc = model.getNumCols();
     
    14561472    FILE * fp = fopen(fn.c_str(),"r");
    14571473    if (!fp) {
    1458       // Try in Samples
    1459       fn = "Samples/input.130";
     1474      // Try in Data/Sample
     1475      fn = "Data/Sample/input.130";
    14601476      fp = fopen(fn.c_str(),"r");
    14611477    }
    14621478    if (!fp) {
    1463       fprintf(stderr,"Unable to open file input.130 in mpsDir or Samples directory\n");
     1479      fprintf(stderr,"Unable to open file input.130 in mpsDir or Data/Sample directory\n");
    14641480    } else {
    14651481      int problem;
Note: See TracChangeset for help on using the changeset viewer.