Changeset 1376 for trunk/Clp


Ignore:
Timestamp:
Jul 2, 2009 12:26:49 PM (10 years ago)
Author:
forrest
Message:

chnages for speed and alternate factorization

Location:
trunk/Clp/src
Files:
12 edited

Legend:

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

    r1370 r1376  
    672672                                        double & objectiveChange)
    673673{
    674   double * work = primalUpdate->denseVector();
     674  double * COIN_RESTRICT work = primalUpdate->denseVector();
    675675  int number = primalUpdate->getNumElements();
    676   int * which = primalUpdate->getIndices();
     676  int * COIN_RESTRICT which = primalUpdate->getIndices();
    677677  int i;
    678678  double changeObj=0.0;
    679679  double tolerance=model_->currentPrimalTolerance();
    680   const int * pivotVariable = model_->pivotVariable();
    681   double * infeas = infeasible_->denseVector();
    682   int pivotRow = model_->pivotRow();
    683   double * solution = model_->solutionRegion();
     680  const int * COIN_RESTRICT pivotVariable = model_->pivotVariable();
     681  double * COIN_RESTRICT infeas = infeasible_->denseVector();
     682  double * COIN_RESTRICT solution = model_->solutionRegion();
     683  const double * COIN_RESTRICT costModel = model_->costRegion();
     684  const double * COIN_RESTRICT lowerModel = model_->lowerRegion();
     685  const double * COIN_RESTRICT upperModel = model_->upperRegion();
    684686#ifdef COLUMN_BIAS
    685687  int numberColumns = model_->numberColumns();
     
    690692      int iPivot=pivotVariable[iRow];
    691693      double value = solution[iPivot];
    692       double cost = model_->cost(iPivot);
     694      double cost = costModel[iPivot];
    693695      double change = primalRatio*work[i];
    694696      work[i]=0.0;
    695697      value -= change;
    696698      changeObj -= change*cost;
     699      double lower = lowerModel[iPivot];
     700      double upper = upperModel[iPivot];
    697701      solution[iPivot] = value;
    698       double lower = model_->lower(iPivot);
    699       double upper = model_->upper(iPivot);
    700       // But if pivot row then use value of incoming
    701       // Although it is safer to recompute before next selection
    702       // in case something odd happens
    703       if (iRow==pivotRow) {
    704         iPivot = model_->sequenceIn();
    705         lower = model_->lower(iPivot);
    706         upper = model_->upper(iPivot);
    707         value = model_->valueIncomingDual();
    708       }
    709702      if (value<lower-tolerance) {
    710703        value -= lower;
     
    750743      int iPivot=pivotVariable[iRow];
    751744      double value = solution[iPivot];
    752       double cost = model_->cost(iPivot);
     745      double cost = costModel[iPivot];
    753746      double change = primalRatio*work[iRow];
    754747      value -= change;
    755748      changeObj -= change*cost;
     749      double lower = lowerModel[iPivot];
     750      double upper = upperModel[iPivot];
    756751      solution[iPivot] = value;
    757       double lower = model_->lower(iPivot);
    758       double upper = model_->upper(iPivot);
    759       // But if pivot row then use value of incoming
    760       // Although it is safer to recompute before next selection
    761       // in case something odd happens
    762       if (iRow==pivotRow) {
    763         iPivot = model_->sequenceIn();
    764         lower = model_->lower(iPivot);
    765         upper = model_->upper(iPivot);
    766         value = model_->valueIncomingDual();
    767       }
    768752      if (value<lower-tolerance) {
    769753        value -= lower;
     
    805789      work[iRow]=0.0;
    806790    }
     791  }
     792  // Do pivot row
     793  {
     794    int iRow = model_->pivotRow();
     795    // feasible - was it infeasible - if so set tiny
     796    assert (infeas[iRow]);
     797    //if (infeas[iRow])
     798    infeas[iRow] = COIN_INDEXED_REALLY_TINY_ELEMENT;
    807799  }
    808800  primalUpdate->setNumElements(0);
  • trunk/Clp/src/ClpFactorization.cpp

    r1371 r1376  
    723723                                  int pivotRow,
    724724                                  double pivotCheck ,
    725                                   bool checkBeforeModifying)
     725                                  bool checkBeforeModifying,
     726                                  double acceptablePivot)
    726727{
    727728  int returnCode;
     
    737738                                              pivotRow,
    738739                                              pivotCheck,
    739                                               checkBeforeModifying);
     740                                                   checkBeforeModifying,
     741                                                   acceptablePivot);
    740742    } else {
    741743      returnCode= CoinFactorization::replaceColumnPFI(tableauColumn,
     
    752754                                                        pivotRow,
    753755                                                        pivotCheck,
    754                                                         checkBeforeModifying);
     756                                                    checkBeforeModifying,
     757                                                    acceptablePivot);
    755758      networkBasis_->replaceColumn(regionSparse,
    756759                                   pivotRow);
     
    11081111#else
    11091112// This one allows multiple factorizations
     1113#if CLP_MULTIPLE_FACTORIZATIONS == 1
     1114typedef CoinDenseFactorization CoinOtherFactorization;
     1115typedef CoinOslFactorization CoinOtherFactorization;
     1116#elif CLP_MULTIPLE_FACTORIZATIONS == 2
     1117#include "CoinSimpFactorization.hpp"
     1118typedef CoinSimpFactorization CoinOtherFactorization;
     1119typedef CoinOslFactorization CoinOtherFactorization;
     1120#elif CLP_MULTIPLE_FACTORIZATIONS == 3
     1121#include "CoinSimpFactorization.hpp"
     1122#define CoinOslFactorization CoinDenseFactorization
     1123#elif CLP_MULTIPLE_FACTORIZATIONS == 4
     1124#include "CoinDenseFactorization.hpp"
     1125#include "CoinSimpFactorization.hpp"
     1126#include "CoinOslFactorization.hpp"
     1127#endif
    11101128
    11111129//-------------------------------------------------------------------
     
    11201138  coinFactorizationA_ = new CoinFactorization() ;
    11211139  coinFactorizationB_ = NULL;
    1122   //coinFactorizationB_ = new CoinSmallFactorization();
     1140  //coinFactorizationB_ = new CoinOtherFactorization();
    11231141  forceB_=0;
    11241142  goOslThreshold_ = -1;
     
    11471165  goSmallThreshold_ = rhs.goSmallThreshold_;
    11481166  int goDense = 0;
    1149   if (denseIfSmaller>0&&!rhs.coinFactorizationB_) {
     1167  if (denseIfSmaller>0&&denseIfSmaller<=goDenseThreshold_) {
     1168    CoinDenseFactorization * denseR =
     1169      dynamic_cast<CoinDenseFactorization *>(rhs.coinFactorizationB_);
     1170    if (!denseR)
     1171      goDense=1;
     1172  }
     1173  if (denseIfSmaller>0&&!rhs.coinFactorizationB_) {
    11501174    if (denseIfSmaller<=goDenseThreshold_)
    11511175      goDense=1;
     
    12151239}
    12161240
    1217 ClpFactorization::ClpFactorization (const CoinSmallFactorization & rhs)
     1241ClpFactorization::ClpFactorization (const CoinOtherFactorization & rhs)
    12181242{
    12191243#ifdef CLP_FACTORIZATION_INSTRUMENT
     
    12251249  coinFactorizationA_ = NULL;
    12261250  coinFactorizationB_ = rhs.clone();
    1227   //coinFactorizationB_ = new CoinSmallFactorization(rhs);
     1251  //coinFactorizationB_ = new CoinOtherFactorization(rhs);
    12281252  forceB_=0;
    12291253  goOslThreshold_ = -1;
     
    19651989        //int * pivotTemp = pivotColumn_.array();
    19661990        //ClpDisjointCopyN ( pivotVariable, numberRows , pivotTemp  );
     1991#ifndef NDEBUG
     1992        CoinFillN(pivotVariable,numberRows,-1);
     1993#endif
    19671994        // Redo pivot order
    19681995        for (i=0;i<numberRowBasic;i++) {
    19691996          int k = pivotTemp[i];
    19701997          // so rowIsBasic[k] would be permuteBack[back[i]]
    1971           pivotVariable[permuteBack[back[i]]]=k+numberColumns;
     1998          int j=permuteBack[back[i]];
     1999          assert (pivotVariable[j]==-1);
     2000          pivotVariable[j]=k+numberColumns;
    19722001        }
    19732002        for (;i<useNumberRows;i++) {
    19742003          int k = pivotTemp[i];
    19752004          // so rowIsBasic[k] would be permuteBack[back[i]]
    1976           pivotVariable[permuteBack[back[i]]]=k;
     2005          int j=permuteBack[back[i]];
     2006          assert (pivotVariable[j]==-1);
     2007          pivotVariable[j]=k;
    19772008        }
    19782009#if 0
     
    22282259                                  int pivotRow,
    22292260                                  double pivotCheck ,
    2230                                   bool checkBeforeModifying)
     2261                                  bool checkBeforeModifying,
     2262                                       double acceptablePivot)
    22312263{
    22322264#ifndef SLIM_CLP
     
    22432275                                                        pivotRow,
    22442276                                                        pivotCheck,
    2245                                                         checkBeforeModifying);
     2277                                                        checkBeforeModifying,
     2278                                                        acceptablePivot);
    22462279      } else {
    22472280        bool tab = coinFactorizationB_->wantsTableauColumn();
     
    22522285                                             pivotRow,
    22532286                                             pivotCheck,
    2254                                              checkBeforeModifying);
     2287                                             checkBeforeModifying,
     2288                                             acceptablePivot);
    22552289#ifdef CLP_DEBUG
    22562290        // check basic
     
    22922326                                                pivotRow,
    22932327                                                pivotCheck,
    2294                                                 checkBeforeModifying);
     2328                                                          checkBeforeModifying,
     2329                                                          acceptablePivot);
    22952330      networkBasis_->replaceColumn(regionSparse,
    22962331                                   pivotRow);
  • trunk/Clp/src/ClpFactorization.hpp

    r1371 r1376  
    1212class ClpSimplex;
    1313class ClpNetworkBasis;
     14class CoinOtherFactorization;
    1415#ifndef CLP_MULTIPLE_FACTORIZATIONS
    1516#ifdef CLP_OSL
     
    1920#endif
    2021#endif   
    21 #if CLP_MULTIPLE_FACTORIZATIONS == 1
     22#ifdef CLP_MULTIPLE_FACTORIZATIONS
    2223#include "CoinDenseFactorization.hpp"
    23 typedef CoinDenseFactorization CoinSmallFactorization;
    24 typedef CoinOslFactorization CoinSmallFactorization;
    25 #elif CLP_MULTIPLE_FACTORIZATIONS == 2
    26 #include "CoinSimpFactorization.hpp"
    27 typedef CoinSimpFactorization CoinSmallFactorization;
    28 typedef CoinOslFactorization CoinSmallFactorization;
    29 #elif CLP_MULTIPLE_FACTORIZATIONS == 3
    30 #include "CoinDenseFactorization.hpp"
    31 #include "CoinSimpFactorization.hpp"
    32 #define CoinOslFactorization CoinDenseFactorization
    33 #elif CLP_MULTIPLE_FACTORIZATIONS == 4
    34 #include "CoinDenseFactorization.hpp"
    35 #include "CoinSimpFactorization.hpp"
    36 #include "CoinOslFactorization.hpp"
    3724#endif
    3825
     
    8067  ClpFactorization(const ClpFactorization&,int denseIfSmaller=0);
    8168#ifdef CLP_MULTIPLE_FACTORIZATIONS   
    82    /** The copy constructor from an CoinSmallFactorization. */
    83    ClpFactorization(const CoinSmallFactorization&);
     69   /** The copy constructor from an CoinOtherFactorization. */
     70   ClpFactorization(const CoinOtherFactorization&);
    8471#endif
    8572   ClpFactorization& operator=(const ClpFactorization&);
     
    10289                      int pivotRow,
    10390                      double pivotCheck ,
    104                       bool checkBeforeModifying=false);
     91                      bool checkBeforeModifying=false,
     92                      double acceptablePivot=1.0e-8);
    10593  //@}
    10694
     
    226214  /// Level of detail of messages
    227215  inline int messageLevel (  ) const {
    228     if (coinFactorizationA_) return coinFactorizationA_->messageLevel();  else return 0 ;
     216    if (coinFactorizationA_) return coinFactorizationA_->messageLevel();  else return 1 ;
    229217  }
    230218  /// Set level of detail of messages
     
    234222  /// Get rid of all memory
    235223  inline void clearArrays()
    236   { if (coinFactorizationA_) coinFactorizationA_->clearArrays();}
     224  { if (coinFactorizationA_)
     225      coinFactorizationA_->clearArrays();
     226    else if (coinFactorizationB_)
     227      coinFactorizationB_->clearArrays();
     228  }
    237229  /// Number of Rows after factorization
    238230  inline int numberRows (  ) const {
     
    247239  /// Pivot tolerance
    248240  inline double pivotTolerance (  ) const {
    249     if (coinFactorizationA_) return coinFactorizationA_->pivotTolerance();  else return 1.0e-8 ;
     241    if (coinFactorizationA_) return coinFactorizationA_->pivotTolerance();  else if (coinFactorizationB_) return coinFactorizationB_->pivotTolerance();return 1.0e-8 ;
    250242  }
    251243  /// Set pivot tolerance
    252244  inline void pivotTolerance (  double value) {
    253245    if (coinFactorizationA_) coinFactorizationA_->pivotTolerance(value);
     246    else if (coinFactorizationB_) coinFactorizationB_->pivotTolerance(value);
    254247  }
    255248  /// Allows change of pivot accuracy check 1.0 == none >1.0 relaxed
     
    267260  /// Delete all stuff (leaves as after CoinFactorization())
    268261  inline void almostDestructor()
    269   { if (coinFactorizationA_) coinFactorizationA_->almostDestructor(); }
     262  { if (coinFactorizationA_)
     263      coinFactorizationA_->almostDestructor();
     264    else if (coinFactorizationB_)
     265      coinFactorizationB_->clearArrays();
     266  }
    270267  /// Returns areaFactor but adjusted for dense
    271268  inline double adjustedAreaFactor() const
     
    366363  /// Pointer to CoinFactorization
    367364  CoinFactorization * coinFactorizationA_;
    368   /// Pointer to CoinSmallFactorization
    369   CoinSmallFactorization * coinFactorizationB_;
     365  /// Pointer to CoinOtherFactorization
     366  CoinOtherFactorization * coinFactorizationB_;
    370367  /// If nonzero force use of 1,dense 2,small 3,osl
    371368  int forceB_;
  • trunk/Clp/src/ClpHelperFunctions.hpp

    r1370 r1376  
    3737{
    3838  for (int i=0;i<size;i++)
    39     to[i]=from[i];
     39    to[i]=static_cast<double>(from[i]);
    4040}
    4141inline CoinWorkDouble
  • trunk/Clp/src/ClpModel.hpp

    r1370 r1376  
    850850      32768 - called from Osi
    851851      65536 - keep arrays around as much as possible (also use maximumR/C)
    852       131072 - unused
     852      131072 - transposeTimes is -1.0 and can skip basic and fixed
    853853      262144 - extra copy of scaled matrix
    854854      524288 - Clp fast dual
  • trunk/Clp/src/ClpPackedMatrix.cpp

    r1370 r1376  
    849849      }
    850850      if (!columnCopy_) {
    851         numberNonZero=thisMatrix->gutsOfTransposeTimesUnscaled(pi,
    852                                                                columnArray->getIndices(),
    853                                                                columnArray->denseVector(),
    854                                                                zeroTolerance);
     851        if ((model->specialOptions(),131072)!=0) {
     852          if(model->spareIntArray_[0]>0) {
     853            CoinIndexedVector * spareArray = model->rowArray(3);
     854            // also do dualColumn stuff
     855            double * spare = spareArray->denseVector();
     856            int * spareIndex = spareArray->getIndices();
     857            const double * reducedCost=model->djRegion(0);
     858            double multiplier[]={-1.0,1.0};
     859            double dualT = - model->currentDualTolerance();
     860            double acceptablePivot=model->spareDoubleArray_[0];
     861            // We can also see if infeasible or pivoting on free
     862            double tentativeTheta = 1.0e15;
     863            double upperTheta = 1.0e31;
     864            double bestPossible=0.0;
     865            int addSequence=model->numberColumns();
     866            const unsigned char * statusArray = model->statusArray()+addSequence;
     867            int numberRemaining=0;
     868            assert (scalar==-1.0);
     869            for (i=0;i<numberInRowArray;i++) {
     870              int iSequence = whichRow[i];
     871              int iStatus = (statusArray[iSequence]&3)-1;
     872              if (iStatus) {
     873                double mult=multiplier[iStatus-1];
     874                double alpha=piOld[i]*mult;
     875                double oldValue;
     876                double value;
     877                if (alpha>0.0) {
     878                  oldValue = reducedCost[iSequence]*mult;
     879                  value = oldValue-tentativeTheta*alpha;
     880                  if (value<dualT) {
     881                    bestPossible = CoinMax(bestPossible,alpha);
     882                    value = oldValue-upperTheta*alpha;
     883                    if (value<dualT && alpha>=acceptablePivot) {
     884                      upperTheta = (oldValue-dualT)/alpha;
     885                      //tentativeTheta = CoinMin(2.0*upperTheta,tentativeTheta);
     886                    }
     887                    // add to list
     888                    spare[numberRemaining]=alpha*mult;
     889                    spareIndex[numberRemaining++]=iSequence+addSequence;
     890                  }
     891                }
     892              }
     893            }
     894            numberNonZero=
     895              thisMatrix->gutsOfTransposeTimesUnscaled(pi,
     896                                                       columnArray->getIndices(),
     897                                                       columnArray->denseVector(),
     898                                                       model->statusArray(),
     899                                                       spareIndex,
     900                                                       spare,
     901                                                       model->djRegion(1),
     902                                                       upperTheta,
     903                                                       bestPossible,
     904                                                       acceptablePivot,
     905                                                       model->currentDualTolerance(),
     906                                                       numberRemaining,
     907                                                       zeroTolerance);
     908            model->spareDoubleArray_[0]=upperTheta;
     909            model->spareDoubleArray_[1]=bestPossible;
     910            spareArray->setNumElements(numberRemaining);
     911            // signal partially done
     912            model->spareIntArray_[0]=-2;
     913          } else {
     914            numberNonZero=
     915              thisMatrix->gutsOfTransposeTimesUnscaled(pi,
     916                                                       columnArray->getIndices(),
     917                                                       columnArray->denseVector(),
     918                                                       model->statusArray(),
     919                                                       zeroTolerance);
     920          }
     921        } else {
     922          numberNonZero=
     923            thisMatrix->gutsOfTransposeTimesUnscaled(pi,
     924                                                   columnArray->getIndices(),
     925                                                   columnArray->denseVector(),
     926                                                   zeroTolerance);
     927        }
    855928        columnArray->setNumElements(numberNonZero);
    856929        //xA++;
     
    881954      const double * columnScale = model->columnScale();
    882955      if (!columnCopy_) {
    883         numberNonZero=gutsOfTransposeTimesScaled(pi,columnScale,
    884                                                  columnArray->getIndices(),
    885                                                  columnArray->denseVector(),
    886                                                  zeroTolerance);
     956        if ((model->specialOptions(),131072)!=0)
     957          numberNonZero=
     958            gutsOfTransposeTimesScaled(pi,columnScale,
     959                                       columnArray->getIndices(),
     960                                       columnArray->denseVector(),
     961                                       model->statusArray(),
     962                                       zeroTolerance);
     963        else
     964          numberNonZero=
     965            gutsOfTransposeTimesScaled(pi,columnScale,
     966                                       columnArray->getIndices(),
     967                                       columnArray->denseVector(),
     968                                       zeroTolerance);
    887969        columnArray->setNumElements(numberNonZero);
    888970        //xC++;
     
    12911373  return numberNonZero;
    12921374}
     1375// Meat of transposeTimes by column when not scaled
     1376int
     1377ClpPackedMatrix::gutsOfTransposeTimesUnscaled(const double * COIN_RESTRICT pi,
     1378                                    int * COIN_RESTRICT index,
     1379                                    double * COIN_RESTRICT array,
     1380                                              const unsigned char * COIN_RESTRICT status,                     
     1381                                    const double zeroTolerance) const
     1382{
     1383  int numberNonZero=0;
     1384  // get matrix data pointers
     1385  const int * COIN_RESTRICT row = matrix_->getIndices();
     1386  const CoinBigIndex * COIN_RESTRICT columnStart = matrix_->getVectorStarts();
     1387  const double * COIN_RESTRICT elementByColumn = matrix_->getElements();
     1388  double value = 0.0;
     1389  int jColumn=-1;
     1390  for (int iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
     1391    bool wanted = ((status[iColumn]&3)!=1);
     1392    if (fabs(value)>zeroTolerance) {
     1393      array[numberNonZero]=value;
     1394      index[numberNonZero++]=jColumn;
     1395    }
     1396    value = 0.0;
     1397    if (wanted) {
     1398      CoinBigIndex start = columnStart[iColumn];
     1399      CoinBigIndex end = columnStart[iColumn+1];
     1400      jColumn=iColumn;
     1401      int n=end-start;
     1402      bool odd = (n&1)!=0;
     1403      n = n>>1;
     1404      const int * COIN_RESTRICT rowThis = row+start;
     1405      const double * COIN_RESTRICT elementThis = elementByColumn+start;
     1406      for (;n;n--) {
     1407        int iRow0 = *rowThis;
     1408        int iRow1 = *(rowThis+1);
     1409        rowThis+=2;
     1410        value += pi[iRow0]*(*elementThis);
     1411        value += pi[iRow1]*(*(elementThis+1));
     1412        elementThis+=2;
     1413      }
     1414      if (odd) {
     1415        int iRow = *rowThis;
     1416        value += pi[iRow]*(*elementThis);
     1417      }
     1418    }
     1419  }
     1420  if (fabs(value)>zeroTolerance) {
     1421    array[numberNonZero]=value;
     1422    index[numberNonZero++]=jColumn;
     1423  }
     1424  return numberNonZero;
     1425}
     1426/* Meat of transposeTimes by column when not scaled and skipping
     1427   and doing part of dualColumn */
     1428 int
     1429ClpPackedMatrix::gutsOfTransposeTimesUnscaled(const double * COIN_RESTRICT pi,
     1430                                    int * COIN_RESTRICT index,
     1431                                    double * COIN_RESTRICT array,
     1432                                   const unsigned char * COIN_RESTRICT status,
     1433                                    int * COIN_RESTRICT spareIndex,
     1434                                    double * COIN_RESTRICT spareArray,
     1435                                    const double * COIN_RESTRICT reducedCost,
     1436                                   double & upperThetaP,
     1437                                   double & bestPossibleP,
     1438                                   double acceptablePivot,
     1439                                   double dualTolerance,
     1440                                   int & numberRemainingP,
     1441                                    const double zeroTolerance) const
     1442 {
     1443  double tentativeTheta = 1.0e15;
     1444  int numberRemaining = numberRemainingP;
     1445  double upperTheta=upperThetaP;
     1446  double bestPossible=bestPossibleP;
     1447  int numberNonZero=0;
     1448  // get matrix data pointers
     1449  const int * COIN_RESTRICT row = matrix_->getIndices();
     1450  const CoinBigIndex * COIN_RESTRICT columnStart = matrix_->getVectorStarts();
     1451  const double * COIN_RESTRICT elementByColumn = matrix_->getElements();
     1452  double multiplier[]={-1.0,1.0};
     1453  double dualT = - dualTolerance;
     1454  for (int iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
     1455    int wanted = (status[iColumn]&3)-1;
     1456    if (wanted) {
     1457      double value = 0.0;
     1458      CoinBigIndex start = columnStart[iColumn];
     1459      CoinBigIndex end = columnStart[iColumn+1];
     1460      int n=end-start;
     1461      bool odd = (n&1)!=0;
     1462      n = n>>1;
     1463      const int * COIN_RESTRICT rowThis = row+start;
     1464      const double * COIN_RESTRICT elementThis = elementByColumn+start;
     1465      for (;n;n--) {
     1466        int iRow0 = *rowThis;
     1467        int iRow1 = *(rowThis+1);
     1468        rowThis+=2;
     1469        value += pi[iRow0]*(*elementThis);
     1470        value += pi[iRow1]*(*(elementThis+1));
     1471        elementThis+=2;
     1472      }
     1473      if (odd) {
     1474        int iRow = *rowThis;
     1475        value += pi[iRow]*(*elementThis);
     1476      }
     1477      if (fabs(value)>zeroTolerance) {
     1478        double mult=multiplier[wanted-1];
     1479        double alpha = value*mult;
     1480        array[numberNonZero]=value;
     1481        index[numberNonZero++]=iColumn;
     1482        if (alpha>0.0) {
     1483          double oldValue = reducedCost[iColumn]*mult;
     1484          double value = oldValue-tentativeTheta*alpha;
     1485          if (value<dualT) {
     1486            bestPossible = CoinMax(bestPossible,alpha);
     1487            value = oldValue-upperTheta*alpha;
     1488            if (value<dualT && alpha>=acceptablePivot) {
     1489              upperTheta = (oldValue-dualT)/alpha;
     1490              //tentativeTheta = CoinMin(2.0*upperTheta,tentativeTheta);
     1491            }
     1492            // add to list
     1493            spareArray[numberRemaining]=alpha*mult;
     1494            spareIndex[numberRemaining++]=iColumn;
     1495          }
     1496        }
     1497      }
     1498    }
     1499  }
     1500  numberRemainingP = numberRemaining;
     1501  upperThetaP=upperTheta;
     1502  bestPossibleP=bestPossible;
     1503  return numberNonZero;
     1504 }
     1505// Meat of transposeTimes by column when scaled
     1506int
     1507ClpPackedMatrix::gutsOfTransposeTimesScaled(const double * COIN_RESTRICT pi,
     1508                                 const double * COIN_RESTRICT columnScale,
     1509                                 int * COIN_RESTRICT index,
     1510                                 double * COIN_RESTRICT array,
     1511                                            const unsigned char * COIN_RESTRICT status,                          const double zeroTolerance) const
     1512{
     1513  int numberNonZero=0;
     1514  // get matrix data pointers
     1515  const int * COIN_RESTRICT row = matrix_->getIndices();
     1516  const CoinBigIndex * COIN_RESTRICT columnStart = matrix_->getVectorStarts();
     1517  const double * COIN_RESTRICT elementByColumn = matrix_->getElements();
     1518  double value = 0.0;
     1519  int jColumn=-1;
     1520  for (int iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
     1521    bool wanted = ((status[iColumn]&3)!=1);
     1522    if (fabs(value)>zeroTolerance) {
     1523      array[numberNonZero]=value;
     1524      index[numberNonZero++]=jColumn;
     1525    }
     1526    value = 0.0;
     1527    if (wanted) {
     1528      double scale = columnScale[iColumn];
     1529      CoinBigIndex start = columnStart[iColumn];
     1530      CoinBigIndex end = columnStart[iColumn+1];
     1531      jColumn=iColumn;
     1532      for (CoinBigIndex j=start; j<end;j++) {
     1533        int iRow = row[j];
     1534        value += pi[iRow]*elementByColumn[j];
     1535      }
     1536      value *= scale;
     1537    }
     1538  }
     1539  if (fabs(value)>zeroTolerance) {
     1540    array[numberNonZero]=value;
     1541    index[numberNonZero++]=jColumn;
     1542  }
     1543  return numberNonZero;
     1544}
    12931545// Meat of transposeTimes by row n > K if packed - returns number nonzero
    12941546int
     
    13101562    int iRow = whichRow[i];
    13111563    double value = pi[i]*scalar;
    1312     CoinBigIndex j;
     1564    CoinBigIndex start=rowStart[iRow];
     1565    CoinBigIndex end=rowStart[iRow+1];
     1566    int n=end-start;
     1567    const int * COIN_RESTRICT columnThis = column+start;
     1568    const double * COIN_RESTRICT elementThis = element+start;
     1569
    13131570    // could do by twos
    1314     for (j=rowStart[iRow];j<rowStart[iRow+1];j++) {
    1315       int iColumn = column[j];
    1316       double elValue = element[j];
     1571    for (;n;n--) {
     1572      int iColumn = *columnThis;
     1573      columnThis++;
     1574      double elValue = *elementThis;
     1575      elementThis++;
    13171576      elValue *= value;
    13181577      output[iColumn] += elValue;
     
    28093068}
    28103069#endif
     3070//#define SQRT_ARRAY
     3071#ifdef SQRT_ARRAY
     3072static void doSqrts(double * array,int n)
     3073{
     3074  for (int i=0;i<n;i++)
     3075    array[i] = 1.0/sqrt(array[i]);
     3076}
     3077#endif
    28113078//static int scale_stats[5]={0,0,0,0,0};
    28123079// Creates scales for column copy (rowCopy in model may be modified)
     
    30393306          for (iRow=0;iRow<numberRows;iRow++) {
    30403307            CoinBigIndex j;
    3041             largest=1.0e-20;
     3308            largest=1.0e-50;
    30423309            smallest=1.0e50;
    30433310            for (j=rowStart[iRow];j<rowStart[iRow+1];j++) {
     
    30513318            }
    30523319           
     3320#ifdef SQRT_ARRAY
     3321            rowScale[iRow]=smallest*largest;
     3322#else
    30533323            rowScale[iRow]=1.0/sqrt(smallest*largest);
     3324#endif
    30543325            //rowScale[iRow]=CoinMax(1.0e-10,CoinMin(1.0e10,rowScale[iRow]));
    30553326            if (extraDetails) {
     
    30583329            }
    30593330          }
     3331#ifdef SQRT_ARRAY
     3332          doSqrts(rowScale,numberRows);
     3333#endif
    30603334#ifdef USE_OBJECTIVE
    30613335          largest=1.0e-20;
     
    30823356            if (usefulColumn[iColumn]) {
    30833357              CoinBigIndex j;
    3084               largest=1.0e-20;
     3358              largest=1.0e-50;
    30853359              smallest=1.0e50;
    30863360              for (j=columnStart[iColumn];
     
    30993373              }
    31003374#endif
     3375#ifdef SQRT_ARRAY
     3376              columnScale[iColumn]=smallest*largest;
     3377#else
    31013378              columnScale[iColumn]=1.0/sqrt(smallest*largest);
     3379#endif
    31023380              //columnScale[iColumn]=CoinMax(1.0e-10,CoinMin(1.0e10,columnScale[iColumn]));
    31033381            }
    31043382          }
     3383#ifdef SQRT_ARRAY
     3384          doSqrts(columnScale,numberColumns);
     3385#endif
    31053386        }
    31063387#ifdef USE_OBJECTIVE
  • trunk/Clp/src/ClpPackedMatrix.hpp

    r1370 r1376  
    358358                                 double * COIN_RESTRICT array,
    359359                                 const double tolerance) const;
     360  /// Meat of transposeTimes by column when not scaled and skipping
     361  int gutsOfTransposeTimesUnscaled(const double * COIN_RESTRICT pi,
     362                                    int * COIN_RESTRICT index,
     363                                    double * COIN_RESTRICT array,
     364                                   const unsigned char * status,
     365                                    const double tolerance) const;
     366  /** Meat of transposeTimes by column when not scaled and skipping
     367      and doing part of dualColumn */
     368  int gutsOfTransposeTimesUnscaled(const double * COIN_RESTRICT pi,
     369                                    int * COIN_RESTRICT index,
     370                                    double * COIN_RESTRICT array,
     371                                   const unsigned char * status,
     372                                    int * COIN_RESTRICT spareIndex,
     373                                    double * COIN_RESTRICT spareArray,
     374                                    const double * COIN_RESTRICT reducedCost,
     375                                   double & upperTheta,
     376                                   double & bestPossible,
     377                                   double acceptablePivot,
     378                                   double dualTolerance,
     379                                   int & numberRemaining,
     380                                    const double zeroTolerance) const;
     381  /// Meat of transposeTimes by column when scaled and skipping
     382  int gutsOfTransposeTimesScaled(const double * COIN_RESTRICT pi,
     383                                 const double * COIN_RESTRICT columnScale,
     384                                 int * COIN_RESTRICT index,
     385                                 double * COIN_RESTRICT array,
     386                                   const unsigned char * status,
     387                                 const double tolerance) const;
    360388  /// Meat of transposeTimes by row n > K if packed - returns number nonzero
    361389  int gutsOfTransposeTimesByRowGEK(const CoinIndexedVector * COIN_RESTRICT piVector,
    362390                                   int * COIN_RESTRICT index,
    363                                    double * COIN_RESTRICT output,
     391                                   double * COIN_RESTRICT output,
    364392                                   int numberColumns,
    365393                                   const double tolerance,
  • trunk/Clp/src/ClpSimplex.cpp

    r1371 r1376  
    4343  columnPrimalSequence_(-2),
    4444  rowPrimalSequence_(-2),
    45   columnDualInfeasibility_(0.0),
    46   rowDualInfeasibility_(0.0),
     45  bestObjectiveValue_(-COIN_DBL_MAX),
    4746  moreSpecialOptions_(2),
    4847  baseIteration_(0),
    4948  primalToleranceToGetOptimal_(-1.0),
    50   remainingDualInfeasibility_(0.0),
    5149  largeValue_(1.0e15),
    5250  largestPrimalError_(0.0),
     
    156154    columnPrimalSequence_(-2),
    157155    rowPrimalSequence_(-2),
    158     columnDualInfeasibility_(0.0),
    159     rowDualInfeasibility_(0.0),
     156    bestObjectiveValue_(-COIN_DBL_MAX),
    160157    moreSpecialOptions_(2),
    161158    baseIteration_(0),
    162159    primalToleranceToGetOptimal_(-1.0),
    163     remainingDualInfeasibility_(0.0),
    164160    largeValue_(1.0e15),
    165161    largestPrimalError_(0.0),
     
    306302    columnPrimalSequence_(-2),
    307303    rowPrimalSequence_(-2),
    308     columnDualInfeasibility_(0.0),
    309     rowDualInfeasibility_(0.0),
     304    bestObjectiveValue_(-COIN_DBL_MAX),
    310305    moreSpecialOptions_(2),
    311306    baseIteration_(0),
    312307    primalToleranceToGetOptimal_(-1.0),
    313     remainingDualInfeasibility_(0.0),
    314308    largeValue_(1.0e15),
    315309    largestPrimalError_(0.0),
     
    327321  upperOut_(-1),
    328322  dualOut_(-1),
    329   dualTolerance_(0.0),
    330   primalTolerance_(0.0),
     323  dualTolerance_(rhs->dualTolerance_),
     324  primalTolerance_(rhs->primalTolerance_),
    331325  sumDualInfeasibilities_(0.0),
    332326  sumPrimalInfeasibilities_(0.0),
     
    608602          double difference= save[iRow];
    609603          if (difference>1.0e-4) {
    610             sort[numberOut]=iPivot;
    611             save[numberOut++]=difference;
     604            sort[numberOut]=iRow;
     605            save[numberOut++]=-difference;
    612606            if (getStatus(iPivot)==basic)
    613607              numberBasic++;
     
    618612        //printf("no errors on basic - going to all slack - numberOut %d\n",numberOut);
    619613        allSlackBasis(true);
    620       }
    621       CoinSort_2(save, save + numberOut, sort,
    622                  CoinFirstGreater_2<double, int>());
     614        CoinIotaN(pivotVariable_,numberRows_,numberColumns_);
     615      }
     616      CoinSort_2(save, save + numberOut, sort);
    623617      numberOut = CoinMin(1000,numberOut);
    624618      for (iRow=0;iRow<numberOut;iRow++) {
    625         int iColumn=sort[iRow];
     619        int jRow=sort[iRow];
     620        int iColumn=pivotVariable_[jRow];
    626621        setColumnStatus(iColumn,superBasic);
     622        setRowStatus(jRow,basic);
     623        pivotVariable_[jRow]=jRow+numberColumns_;
    627624        if (fabs(solution_[iColumn])>1.0e10) {
    628625          if (upper_[iColumn]<0.0) {
     
    16781675    }
    16791676  }
     1677#if 0
     1678  printf("h1 %d %d %g %g %g %g",
     1679      directionOut_
     1680      ,directionIn_,theta_
     1681         ,dualOut_,dualIn_,alpha_);
     1682#endif
    16801683  // change of incoming
    16811684  char rowcol[]={'R','C'};
     
    17421745    handler_->message()<<CoinMessageEol;
    17431746  }
     1747#if 0
     1748  if (numberIterations_>10000)
     1749  printf(" it %d %g %c%d %c%d\n"
     1750      ,numberIterations_,objectiveValue()
     1751      ,rowcol[isColumn(sequenceIn_)],sequenceWithin(sequenceIn_)
     1752         ,rowcol[isColumn(sequenceOut_)],sequenceWithin(sequenceOut_));
     1753#endif
    17441754  if (trustedUserPointer_&&trustedUserPointer_->typeStruct==1) {
    17451755    if (algorithm_>0&&integerType_&&!nonLinearCost_->numberInfeasibilities()) {
     
    19121922  columnPrimalSequence_(-2),
    19131923  rowPrimalSequence_(-2),
    1914   columnDualInfeasibility_(0.0),
    1915   rowDualInfeasibility_(0.0),
     1924  bestObjectiveValue_(rhs.bestObjectiveValue_),
    19161925  moreSpecialOptions_(2),
    19171926  baseIteration_(0),
    19181927  primalToleranceToGetOptimal_(-1.0),
    1919   remainingDualInfeasibility_(0.0),
    19201928  largeValue_(1.0e15),
    19211929  largestPrimalError_(0.0),
     
    20182026  columnPrimalSequence_(-2),
    20192027  rowPrimalSequence_(-2),
    2020   columnDualInfeasibility_(0.0),
    2021   rowDualInfeasibility_(0.0),
     2028  bestObjectiveValue_(-COIN_DBL_MAX),
    20222029  moreSpecialOptions_(2),
    20232030  baseIteration_(0),
    20242031  primalToleranceToGetOptimal_(-1.0),
    2025   remainingDualInfeasibility_(0.0),
    20262032  largeValue_(1.0e15),
    20272033  largestPrimalError_(0.0),
     
    22202226  zeroTolerance_ = rhs.zeroTolerance_;
    22212227  rowPrimalSequence_ = rhs.rowPrimalSequence_;
    2222   columnDualInfeasibility_ = rhs.columnDualInfeasibility_;
    2223   rowDualInfeasibility_ = rhs.rowDualInfeasibility_;
     2228  bestObjectiveValue_ = rhs.bestObjectiveValue_;
    22242229  baseIteration_ = rhs.baseIteration_;
    22252230  primalToleranceToGetOptimal_ = rhs.primalToleranceToGetOptimal_;
    2226   remainingDualInfeasibility_ = rhs.remainingDualInfeasibility_;
    22272231  largeValue_ = rhs.largeValue_;
    22282232  largestPrimalError_ = rhs.largestPrimalError_;
     
    25182522              sumDualInfeasibilities_ += value-dualTolerance_;
    25192523              if (value>possTolerance)
    2520                 bestPossibleImprovement_ += distanceUp*value;
     2524                bestPossibleImprovement_ += CoinMin(distanceUp,1.0e10)*value;
    25212525              if (value>relaxedTolerance)
    25222526                sumOfRelaxedDualInfeasibilities_ += value-relaxedTolerance;
     
    25442548            sumDualInfeasibilities_ += value-dualTolerance_;
    25452549            if (value>possTolerance)
    2546               bestPossibleImprovement_+= value*distanceDown;
     2550              bestPossibleImprovement_+= value*CoinMin(distanceDown,1.0e10);
    25472551            if (value>relaxedTolerance)
    25482552              sumOfRelaxedDualInfeasibilities_ += value-relaxedTolerance;
     
    25802584            sumDualInfeasibilities_ += value-dualTolerance_;
    25812585            if (value>possTolerance)
    2582               bestPossibleImprovement_+= value*distanceUp;
     2586              bestPossibleImprovement_+= value*CoinMin(distanceUp,1.0e10);
    25832587            if (value>relaxedTolerance)
    25842588              sumOfRelaxedDualInfeasibilities_ += value-relaxedTolerance;
     
    25962600            sumDualInfeasibilities_ += value-dualTolerance_;
    25972601            if (value>possTolerance)
    2598               bestPossibleImprovement_+= value*distanceDown;
     2602              bestPossibleImprovement_+= value*CoinMin(distanceDown,1.0e10);
    25992603            if (value>relaxedTolerance)
    26002604              sumOfRelaxedDualInfeasibilities_ += value-relaxedTolerance;
     
    26242628  if ((matrix_->skipDualCheck()&&algorithm_>0&&problemStatus_==-2)||
    26252629      matrix_->rhsOffset(this)) {
     2630    // Say may be free or superbasic
     2631    moreSpecialOptions_ &= ~8;
    26262632    // old way
    26272633    checkPrimalSolution(rowActivityWork_,columnActivityWork_);
     
    26302636  }
    26312637  int iSequence;
    2632 
     2638  assert (dualTolerance_>0.0&&dualTolerance_<1.0e10);
     2639  assert (primalTolerance_>0.0&&primalTolerance_<1.0e10);
    26332640  objectiveValue_ = 0.0;
    26342641  sumPrimalInfeasibilities_=0.0;
     
    26642671
    26652672  int numberTotal = numberRows_ + numberColumns_;
     2673  // Say no free or superbasic
     2674  moreSpecialOptions_ |= 8;
    26662675  for (iSequence=0;iSequence<numberTotal;iSequence++) {
    26672676    double value = solution_[iSequence];
     
    27112720        } else {
    27122721          // may be free
     2722          // Say free or superbasic
     2723          moreSpecialOptions_ &= ~8;
    27132724          djValue *= 0.01;
    27142725          if (fabs(djValue)>dualTolerance) {
     
    36713682      /**********************************************************
    36723683      rowArray_[3] is long enough for rows+columns
     3684      rowArray_[1] is long enough for max(rows,columns)
    36733685      *********************************************************/
    36743686      for (iRow=0;iRow<4;iRow++) {
     
    36763688        if (iRow==3||objective_->type()>1)
    36773689          length += numberColumns_;
     3690        else if (iRow==1)
     3691          length = CoinMax(length,numberColumns_);
    36783692        if ((specialOptions_&65536)==0||!rowArray_[iRow]) {
    36793693          delete rowArray_[iRow];
     
    56815695  barrier.primalDual();
    56825696  int barrierStatus=barrier.status();
    5683   double gap = barrier.complementarityGap();
     5697  double gap = static_cast<double>(barrier.complementarityGap());
    56845698  // get which variables are fixed
    56855699  double * saveLower=NULL;
     
    81658179  otherModel.zeroTolerance_ = zeroTolerance_;
    81668180  otherModel.rowPrimalSequence_ = rowPrimalSequence_;
    8167   otherModel.columnDualInfeasibility_ = columnDualInfeasibility_;
     8181  otherModel.bestObjectiveValue_ = bestObjectiveValue_;
    81688182  otherModel.moreSpecialOptions_ = moreSpecialOptions_;
    8169   otherModel.rowDualInfeasibility_ = rowDualInfeasibility_;
    81708183  otherModel.baseIteration_ = baseIteration_;
    81718184  otherModel.primalToleranceToGetOptimal_ = primalToleranceToGetOptimal_;
    8172   otherModel.remainingDualInfeasibility_ = remainingDualInfeasibility_;
    81738185  otherModel.largestPrimalError_ = largestPrimalError_;
    81748186  otherModel.largestDualError_ = largestDualError_;
  • trunk/Clp/src/ClpSimplex.hpp

    r1371 r1376  
    724724          { alphaAccuracy_ = value;}
    725725public:
     726   /// Objective value
     727   //inline double objectiveValue() const {
     728  //return (objectiveValue_-bestPossibleImprovement_)*optimizationDirection_ - dblParam_[ClpObjOffset];
     729  //}
    726730  /// Set disaster handler
    727731  inline void setDisasterHandler(ClpDisasterHandler * handler)
     
    953957      2 bit - if presolved problem infeasible return
    954958      4 bit - keep arrays like upper_ around
     959      8 bit - no free or superBasic variables
    955960  */
    956961  inline void setMoreSpecialOptions(int value)
     
    11761181  /// Sequence of worst (-1 if feasible)
    11771182  int rowPrimalSequence_;
    1178   /// Worst column dual infeasibility
    1179   double columnDualInfeasibility_;
    1180   /// Worst row dual infeasibility
    1181   double rowDualInfeasibility_;
     1183  /// "Best" objective value
     1184  double bestObjectiveValue_;
    11821185  /// More special options - see set for details
    11831186  int moreSpecialOptions_;
     
    11861189  /// Primal tolerance needed to make dual feasible (<largeTolerance)
    11871190  double primalToleranceToGetOptimal_;
    1188   /// Remaining largest dual infeasibility
    1189   double remainingDualInfeasibility_;
    11901191  /// Large bound value (for complementarity etc)
    11911192  double largeValue_;
  • trunk/Clp/src/ClpSimplexDual.cpp

    r1371 r1376  
    399399  int nPivots=9999;
    400400#endif
     401  // Start can skip some things in transposeTimes
     402  specialOptions_ |= 131072;
    401403  int lastCleaned=0; // last time objective or bounds cleaned up
    402404 
     
    551553  }
    552554#endif
     555  // Stop can skip some things in transposeTimes
     556  specialOptions_ &= ~131072;
    553557}
    554558int
    555559ClpSimplexDual::dual(int ifValuesPass,int startFinishOptions)
    556560{
     561  bestObjectiveValue_=-COIN_DBL_MAX;
    557562  algorithm_ = -1;
    558563  moreSpecialOptions_ &= ~16; // clear check replaceColumn accuracy
     
    590595  restoreData(data);
    591596  dontFactorizePivots_ = saveDont;
     597  if (problemStatus_==3)
     598    objectiveValue_ = CoinMax(bestObjectiveValue_,objectiveValue_-bestPossibleImprovement_);
    592599  return problemStatus_;
    593600}
     
    11411148      else if (factorization_->pivots())
    11421149        acceptablePivot=acceptablePivot_; // relax
     1150      // But factorizations complain if <1.0e-8
     1151      //acceptablePivot=CoinMax(acceptablePivot,1.0e-8);
    11431152      double bestPossiblePivot=1.0;
    11441153      // get sign for finding row of tableau
     
    11561165        sequenceIn_=-1;
    11571166        // put row of tableau in rowArray[0] and columnArray[0]
     1167        assert (!rowArray_[1]->getNumElements());
    11581168        if (!scaledMatrix_) {
     1169          if ((moreSpecialOptions_&8)!=0&&!rowScale_)
     1170            spareIntArray_[0]=1;
    11591171          matrix_->transposeTimes(this,-1.0,
    1160                                   rowArray_[0],rowArray_[3],columnArray_[0]);
    1161         } else {
     1172                                  rowArray_[0],rowArray_[1],columnArray_[0]);
     1173        } else {         
    11621174          double * saveR = rowScale_;
    11631175          double * saveC = columnScale_;
    11641176          rowScale_=NULL;
    11651177          columnScale_=NULL;
     1178          if ((moreSpecialOptions_&8)!=0)
     1179            spareIntArray_[0]=1;
    11661180          scaledMatrix_->transposeTimes(this,-1.0,
    1167                                   rowArray_[0],rowArray_[3],columnArray_[0]);
     1181                                  rowArray_[0],rowArray_[1],columnArray_[0]);
    11681182          rowScale_=saveR;
    11691183          columnScale_=saveC;
     
    14271441                                                         pivotRow_,
    14281442                                                         alpha_,
    1429                                                      (moreSpecialOptions_&16)!=0);
     1443                                                         (moreSpecialOptions_&16)!=0,
     1444                                                         acceptablePivot);
    14301445        // if no pivots, bad update but reasonable alpha - take and invert
    14311446        if (updateStatus==2&&
     
    14681483            updateDualsInDual(rowArray_[0],columnArray_[0],rowArray_[1],
    14691484                              0.0,objectiveChange,true);
     1485            rowArray_[1]->clear();
     1486            columnArray_[0]->clear();
    14701487            continue;
    14711488          }
     
    21212138  if (!fullRecompute) {
    21222139    int i;
    2123     double * reducedCost = djRegion(0);
    2124     const double * lower = lowerRegion(0);
    2125     const double * upper = upperRegion(0);
    2126     const double * cost = costRegion(0);
    2127     double * work;
     2140    double * COIN_RESTRICT reducedCost = djRegion(0);
     2141    const double * COIN_RESTRICT lower = lowerRegion(0);
     2142    const double * COIN_RESTRICT upper = upperRegion(0);
     2143    const double * COIN_RESTRICT cost = costRegion(0);
     2144    double * COIN_RESTRICT work;
    21282145    int number;
    2129     int * which;
     2146    int * COIN_RESTRICT which;
     2147    const unsigned char * COIN_RESTRICT statusArray = status_+numberColumns_;
    21302148    assert(rowArray->packedMode());
    21312149    work = rowArray->denseVector();
    21322150    number = rowArray->getNumElements();
    21332151    which = rowArray->getIndices();
     2152    double multiplier[]={-1.0,1.0};
    21342153    for (i=0;i<number;i++) {
    21352154      int iSequence = which[i];
    21362155      double alphaI = work[i];
    21372156      work[i]=0.0;
    2138      
    2139       Status status = getStatus(iSequence+numberColumns_);
    2140       // more likely to be at upper bound ?
    2141       if (status==atUpperBound) {
    2142 
     2157      int iStatus = (statusArray[iSequence]&3)-1;
     2158      if (iStatus) {
    21432159        double value = reducedCost[iSequence]-theta*alphaI;
    21442160        reducedCost[iSequence]=value;
    2145 
    2146         if (value>tolerance) {
    2147           // to lower bound (if swap)
     2161        double mult=multiplier[iStatus-1];
     2162        value *= mult;
     2163        if (value<-tolerance) {
     2164          // flipping bounds
     2165          double movement = mult*(lower[iSequence] - upper[iSequence]);
    21482166          which[numberInfeasibilities++]=iSequence;
    2149           double movement = lower[iSequence]-upper[iSequence];
    21502167          assert (fabs(movement)<1.0e30);
    21512168#ifdef CLP_DEBUG
     
    21542171                   0,iSequence,value,alphaI,movement);
    21552172#endif
    2156           changeObj += movement*cost[iSequence];
    2157           outputArray->quickAdd(iSequence,-movement);
    2158         }
    2159       } else if (status==atLowerBound) {
    2160 
    2161         double value = reducedCost[iSequence]-theta*alphaI;
    2162         reducedCost[iSequence]=value;
    2163 
    2164         if (value<-tolerance) {
    2165           // to upper bound
    2166           which[numberInfeasibilities++]=iSequence;
    2167           double movement = upper[iSequence] - lower[iSequence];
    2168           assert (fabs(movement)<1.0e30);
     2173          changeObj -= movement*cost[iSequence];
     2174          outputArray->quickAdd(iSequence,movement);
     2175        }
     2176      }
     2177    }
     2178    // Do columns
     2179    reducedCost = djRegion(1);
     2180    lower = lowerRegion(1);
     2181    upper = upperRegion(1);
     2182    cost = costRegion(1);
     2183    // set number of infeasibilities in row array
     2184    numberRowInfeasibilities=numberInfeasibilities;
     2185    rowArray->setNumElements(numberInfeasibilities);
     2186    numberInfeasibilities=0;
     2187    work = columnArray->denseVector();
     2188    number = columnArray->getNumElements();
     2189    which = columnArray->getIndices();
     2190    if ((moreSpecialOptions_&8)!=0) {
     2191      const unsigned char * COIN_RESTRICT statusArray = status_;
     2192      for (i=0;i<number;i++) {
     2193        int iSequence = which[i];
     2194        double alphaI = work[i];
     2195        work[i]=0.0;
     2196       
     2197        int iStatus = (statusArray[iSequence]&3)-1;
     2198        if (iStatus) {
     2199          double value = reducedCost[iSequence]-theta*alphaI;
     2200          reducedCost[iSequence]=value;
     2201          double mult=multiplier[iStatus-1];
     2202          value *= mult;
     2203          if (value<-tolerance) {
     2204            // flipping bounds
     2205            double movement = mult*(upper[iSequence] - lower[iSequence]);
     2206            which[numberInfeasibilities++]=iSequence;
     2207            assert (fabs(movement)<1.0e30);
    21692208#ifdef CLP_DEBUG
    2170           if ((handler_->logLevel()&32))
    2171             printf("%d %d, new dj %g, alpha %g, movement %g\n",
    2172                    0,iSequence,value,alphaI,movement);
    2173 #endif
    2174           changeObj += movement*cost[iSequence];
    2175           outputArray->quickAdd(iSequence,-movement);
     2209            if ((handler_->logLevel()&32))
     2210              printf("%d %d, new dj %g, alpha %g, movement %g\n",
     2211                     1,iSequence,value,alphaI,movement);
     2212#endif
     2213            changeObj += movement*cost[iSequence];
     2214            matrix_->add(this,outputArray,iSequence,movement);
     2215          }
     2216        }
     2217      }
     2218    } else {
     2219      for (i=0;i<number;i++) {
     2220        int iSequence = which[i];
     2221        double alphaI = work[i];
     2222        work[i]=0.0;
     2223       
     2224        Status status = getStatus(iSequence);
     2225        if (status==atLowerBound) {
     2226          double value = reducedCost[iSequence]-theta*alphaI;
     2227          reducedCost[iSequence]=value;
     2228          double movement=0.0;
     2229         
     2230          if (value<-tolerance) {
     2231            // to upper bound
     2232            which[numberInfeasibilities++]=iSequence;
     2233            movement = upper[iSequence] - lower[iSequence];
     2234            assert (fabs(movement)<1.0e30);
     2235#ifdef CLP_DEBUG
     2236            if ((handler_->logLevel()&32))
     2237              printf("%d %d, new dj %g, alpha %g, movement %g\n",
     2238                     1,iSequence,value,alphaI,movement);
     2239#endif
     2240            changeObj += movement*cost[iSequence];
     2241            matrix_->add(this,outputArray,iSequence,movement);
     2242          }
     2243        } else if (status==atUpperBound) {
     2244          double value = reducedCost[iSequence]-theta*alphaI;
     2245          reducedCost[iSequence]=value;
     2246          double movement=0.0;
     2247         
     2248          if (value>tolerance) {
     2249            // to lower bound (if swap)
     2250            which[numberInfeasibilities++]=iSequence;
     2251            movement = lower[iSequence]-upper[iSequence];
     2252            assert (fabs(movement)<1.0e30);
     2253#ifdef CLP_DEBUG
     2254            if ((handler_->logLevel()&32))
     2255              printf("%d %d, new dj %g, alpha %g, movement %g\n",
     2256                     1,iSequence,value,alphaI,movement);
     2257#endif
     2258            changeObj += movement*cost[iSequence];
     2259            matrix_->add(this,outputArray,iSequence,movement);
     2260          }
     2261        } else if (status==isFree) {
     2262          double value = reducedCost[iSequence]-theta*alphaI;
     2263          reducedCost[iSequence]=value;
    21762264        }
    21772265      }
    21782266    }
    21792267  } else  {
    2180     double * solution = solutionRegion(0);
    2181     double * reducedCost = djRegion(0);
    2182     const double * lower = lowerRegion(0);
    2183     const double * upper = upperRegion(0);
    2184     const double * cost = costRegion(0);
    2185     int * which;
     2268    double * COIN_RESTRICT solution = solutionRegion(0);
     2269    double * COIN_RESTRICT reducedCost = djRegion(0);
     2270    const double * COIN_RESTRICT lower = lowerRegion(0);
     2271    const double * COIN_RESTRICT upper = upperRegion(0);
     2272    const double * COIN_RESTRICT cost = costRegion(0);
     2273    int * COIN_RESTRICT which;
    21862274    which = rowArray->getIndices();
    21872275    int iSequence;
     
    22452333      }
    22462334    }
    2247   }
    2248 
    2249   // Do columns
    2250   if (!fullRecompute) {
    2251     int i;
    2252     double * reducedCost = djRegion(1);
    2253     const double * lower = lowerRegion(1);
    2254     const double * upper = upperRegion(1);
    2255     const double * cost = costRegion(1);
    2256     double * work;
    2257     int number;
    2258     int * which;
    2259     // set number of infeasibilities in row array
    2260     numberRowInfeasibilities=numberInfeasibilities;
    2261     rowArray->setNumElements(numberInfeasibilities);
    2262     numberInfeasibilities=0;
    2263     work = columnArray->denseVector();
    2264     number = columnArray->getNumElements();
    2265     which = columnArray->getIndices();
    2266    
    2267     for (i=0;i<number;i++) {
    2268       int iSequence = which[i];
    2269       double alphaI = work[i];
    2270       work[i]=0.0;
    2271      
    2272       Status status = getStatus(iSequence);
    2273       if (status==atLowerBound) {
    2274         double value = reducedCost[iSequence]-theta*alphaI;
    2275         reducedCost[iSequence]=value;
    2276         double movement=0.0;
    2277 
    2278         if (value<-tolerance) {
    2279           // to upper bound
    2280           which[numberInfeasibilities++]=iSequence;
    2281           movement = upper[iSequence] - lower[iSequence];
    2282           assert (fabs(movement)<1.0e30);
    2283 #ifdef CLP_DEBUG
    2284           if ((handler_->logLevel()&32))
    2285             printf("%d %d, new dj %g, alpha %g, movement %g\n",
    2286                    1,iSequence,value,alphaI,movement);
    2287 #endif
    2288           changeObj += movement*cost[iSequence];
    2289           matrix_->add(this,outputArray,iSequence,movement);
    2290         }
    2291       } else if (status==atUpperBound) {
    2292         double value = reducedCost[iSequence]-theta*alphaI;
    2293         reducedCost[iSequence]=value;
    2294         double movement=0.0;
    2295 
    2296         if (value>tolerance) {
    2297           // to lower bound (if swap)
    2298           which[numberInfeasibilities++]=iSequence;
    2299           movement = lower[iSequence]-upper[iSequence];
    2300           assert (fabs(movement)<1.0e30);
    2301 #ifdef CLP_DEBUG
    2302           if ((handler_->logLevel()&32))
    2303             printf("%d %d, new dj %g, alpha %g, movement %g\n",
    2304                    1,iSequence,value,alphaI,movement);
    2305 #endif
    2306           changeObj += movement*cost[iSequence];
    2307           matrix_->add(this,outputArray,iSequence,movement);
    2308         }
    2309       } else if (status==isFree) {
    2310         double value = reducedCost[iSequence]-theta*alphaI;
    2311         reducedCost[iSequence]=value;
    2312       }
    2313     }
    2314   } else {
    2315     double * solution = solutionRegion(1);
    2316     double * reducedCost = djRegion(1);
    2317     const double * lower = lowerRegion(1);
    2318     const double * upper = upperRegion(1);
    2319     const double * cost = costRegion(1);
    2320     int * which;
     2335    // Do columns
     2336    solution = solutionRegion(1);
     2337    reducedCost = djRegion(1);
     2338    lower = lowerRegion(1);
     2339    upper = upperRegion(1);
     2340    cost = costRegion(1);
    23212341    // set number of infeasibilities in row array
    23222342    numberRowInfeasibilities=numberInfeasibilities;
     
    23242344    numberInfeasibilities=0;
    23252345    which = columnArray->getIndices();
    2326     int iSequence;
    23272346    for (iSequence=0;iSequence<numberColumns_;iSequence++) {
    23282347      double value = reducedCost[iSequence];
     
    23842403    }
    23852404  }
     2405
    23862406#ifdef CLP_DEBUG
    23872407  if (fullRecompute&&numberFake_&&(handler_->logLevel()&16)!=0)
     
    29642984  const double * reducedCost;
    29652985  // We can also see if infeasible or pivoting on free
    2966   double tentativeTheta = 1.0e25;
     2986  double tentativeTheta = 1.0e15;
    29672987  double upperTheta = 1.0e31;
    29682988  double freePivot = acceptablePivot;
     
    29712991  int i;
    29722992  badFree=0.0;
    2973   for (int iSection=0;iSection<2;iSection++) {
    2974 
    2975     int addSequence;
    2976 
    2977     if (!iSection) {
    2978       work = rowArray->denseVector();
    2979       number = rowArray->getNumElements();
    2980       which = rowArray->getIndices();
    2981       reducedCost = rowReducedCost_;
    2982       addSequence = numberColumns_;
    2983     } else {
    2984       work = columnArray->denseVector();
    2985       number = columnArray->getNumElements();
    2986       which = columnArray->getIndices();
    2987       reducedCost = reducedCostWork_;
    2988       addSequence = 0;
    2989     }
    2990    
    2991     for (i=0;i<number;i++) {
    2992       int iSequence = which[i];
    2993       double alpha;
    2994       double oldValue;
    2995       double value;
    2996       bool keep;
    2997 
    2998       switch(getStatus(iSequence+addSequence)) {
     2993  if ((moreSpecialOptions_&8)!=0) {
     2994    // No free or super basic
     2995    double multiplier[]={-1.0,1.0};
     2996    double dualT = - dualTolerance_;
     2997    for (int iSection=0;iSection<2;iSection++) {
     2998     
     2999      int addSequence;
     3000      unsigned char * statusArray;
     3001      if (!iSection) {
     3002        work = rowArray->denseVector();
     3003        number = rowArray->getNumElements();
     3004        which = rowArray->getIndices();
     3005        reducedCost = rowReducedCost_;
     3006        addSequence = numberColumns_;
     3007        statusArray = status_+numberColumns_;
     3008      } else {
     3009        work = columnArray->denseVector();
     3010        number = columnArray->getNumElements();
     3011        which = columnArray->getIndices();
     3012        reducedCost = reducedCostWork_;
     3013        addSequence = 0;
     3014        statusArray = status_;
     3015      }
     3016     
     3017      for (i=0;i<number;i++) {
     3018        int iSequence = which[i];
     3019        double alpha;
     3020        double oldValue;
     3021        double value;
     3022       
     3023        assert (getStatus(iSequence+addSequence)!=isFree
     3024                &&getStatus(iSequence+addSequence)!=superBasic);
     3025        int iStatus = (statusArray[iSequence]&3)-1;
     3026        if (iStatus) {
     3027          double mult=multiplier[iStatus-1];
     3028          alpha = work[i]*mult;
     3029          if (alpha>0.0) {
     3030            oldValue = reducedCost[iSequence]*mult;
     3031            value = oldValue-tentativeTheta*alpha;
     3032            if (value<dualT) {
     3033              bestPossible = CoinMax(bestPossible,alpha);
     3034              value = oldValue-upperTheta*alpha;
     3035              if (value<dualT && alpha>=acceptablePivot) {
     3036                upperTheta = (oldValue-dualT)/alpha;
     3037                //tentativeTheta = CoinMin(2.0*upperTheta,tentativeTheta);
     3038              }
     3039              // add to list
     3040              spare[numberRemaining]=alpha*mult;
     3041              index[numberRemaining++]=iSequence+addSequence;
     3042            }
     3043          }
     3044        }
     3045      }
     3046    }
     3047  } else {
     3048    // some free or super basic
     3049    for (int iSection=0;iSection<2;iSection++) {
     3050     
     3051      int addSequence;
     3052     
     3053      if (!iSection) {
     3054        work = rowArray->denseVector();
     3055        number = rowArray->getNumElements();
     3056        which = rowArray->getIndices();
     3057        reducedCost = rowReducedCost_;
     3058        addSequence = numberColumns_;
     3059      } else {
     3060        work = columnArray->denseVector();
     3061        number = columnArray->getNumElements();
     3062        which = columnArray->getIndices();
     3063        reducedCost = reducedCostWork_;
     3064        addSequence = 0;
     3065      }
     3066     
     3067      for (i=0;i<number;i++) {
     3068        int iSequence = which[i];
     3069        double alpha;
     3070        double oldValue;
     3071        double value;
     3072        bool keep;
     3073       
     3074        switch(getStatus(iSequence+addSequence)) {
    29993075         
    3000       case basic:
    3001       case ClpSimplex::isFixed:
    3002         break;
    3003       case isFree:
    3004       case superBasic:
    3005         alpha = work[i];
    3006         bestPossible = CoinMax(bestPossible,fabs(alpha));
    3007         oldValue = reducedCost[iSequence];
    3008         // If free has to be very large - should come in via dualRow
    3009         //if (getStatus(iSequence+addSequence)==isFree&&fabs(alpha)<1.0e-3)
    3010         //break;
    3011         if (oldValue>dualTolerance_) {
    3012           keep = true;
    3013         } else if (oldValue<-dualTolerance_) {
    3014           keep = true;
    3015         } else {
    3016           if (fabs(alpha)>CoinMax(10.0*acceptablePivot,1.0e-5)) {
     3076        case basic:
     3077        case ClpSimplex::isFixed:
     3078          break;
     3079        case isFree:
     3080        case superBasic:
     3081          alpha = work[i];
     3082          bestPossible = CoinMax(bestPossible,fabs(alpha));
     3083          oldValue = reducedCost[iSequence];
     3084          // If free has to be very large - should come in via dualRow
     3085          //if (getStatus(iSequence+addSequence)==isFree&&fabs(alpha)<1.0e-3)
     3086          //break;
     3087          if (oldValue>dualTolerance_) {
     3088            keep = true;
     3089          } else if (oldValue<-dualTolerance_) {
    30173090            keep = true;
    30183091          } else {
    3019             keep=false;
    3020             badFree=CoinMax(badFree,fabs(alpha));
    3021           }
    3022         }
    3023         if (keep) {
    3024           // free - choose largest
    3025           if (fabs(alpha)>freePivot) {
    3026             freePivot=fabs(alpha);
    3027             sequenceIn_ = iSequence + addSequence;
    3028             theta_=oldValue/alpha;
    3029             alpha_=alpha;
    3030           }
    3031         }
    3032         break;
    3033       case atUpperBound:
    3034         alpha = work[i];
    3035         oldValue = reducedCost[iSequence];
    3036         value = oldValue-tentativeTheta*alpha;
    3037         //assert (oldValue<=dualTolerance_*1.0001);
    3038         if (value>dualTolerance_) {
    3039           bestPossible = CoinMax(bestPossible,-alpha);
    3040           value = oldValue-upperTheta*alpha;
    3041           if (value>dualTolerance_ && -alpha>=acceptablePivot)
    3042             upperTheta = (oldValue-dualTolerance_)/alpha;
    3043           // add to list
    3044           spare[numberRemaining]=alpha;
    3045           index[numberRemaining++]=iSequence+addSequence;
    3046         }
    3047         break;
    3048       case atLowerBound:
    3049         alpha = work[i];
    3050         oldValue = reducedCost[iSequence];
    3051         value = oldValue-tentativeTheta*alpha;
    3052         //assert (oldValue>=-dualTolerance_*1.0001);
    3053         if (value<-dualTolerance_) {
    3054           bestPossible = CoinMax(bestPossible,alpha);
    3055           value = oldValue-upperTheta*alpha;
    3056           if (value<-dualTolerance_ && alpha>=acceptablePivot)
    3057             upperTheta = (oldValue+dualTolerance_)/alpha;
    3058           // add to list
    3059           spare[numberRemaining]=alpha;
    3060           index[numberRemaining++]=iSequence+addSequence;
    3061         }
    3062         break;
     3092            if (fabs(alpha)>CoinMax(10.0*acceptablePivot,1.0e-5)) {
     3093              keep = true;
     3094            } else {
     3095              keep=false;
     3096              badFree=CoinMax(badFree,fabs(alpha));
     3097            }
     3098          }
     3099          if (keep) {
     3100            // free - choose largest
     3101            if (fabs(alpha)>freePivot) {
     3102              freePivot=fabs(alpha);
     3103              sequenceIn_ = iSequence + addSequence;
     3104              theta_=oldValue/alpha;
     3105              alpha_=alpha;
     3106            }
     3107          }
     3108          break;
     3109        case atUpperBound:
     3110          alpha = work[i];
     3111          oldValue = reducedCost[iSequence];
     3112          value = oldValue-tentativeTheta*alpha;
     3113          //assert (oldValue<=dualTolerance_*1.0001);
     3114          if (value>dualTolerance_) {
     3115            bestPossible = CoinMax(bestPossible,-alpha);
     3116            value = oldValue-upperTheta*alpha;
     3117            if (value>dualTolerance_ && -alpha>=acceptablePivot) {
     3118              upperTheta = (oldValue-dualTolerance_)/alpha;
     3119              //tentativeTheta = CoinMin(2.0*upperTheta,tentativeTheta);
     3120            }
     3121            // add to list
     3122            spare[numberRemaining]=alpha;
     3123            index[numberRemaining++]=iSequence+addSequence;
     3124          }
     3125          break;
     3126        case atLowerBound:
     3127          alpha = work[i];
     3128          oldValue = reducedCost[iSequence];
     3129          value = oldValue-tentativeTheta*alpha;
     3130          //assert (oldValue>=-dualTolerance_*1.0001);
     3131          if (value<-dualTolerance_) {
     3132            bestPossible = CoinMax(bestPossible,alpha);
     3133            value = oldValue-upperTheta*alpha;
     3134            if (value<-dualTolerance_ && alpha>=acceptablePivot) {
     3135              upperTheta = (oldValue+dualTolerance_)/alpha;
     3136              //tentativeTheta = CoinMin(2.0*upperTheta,tentativeTheta);
     3137            }
     3138            // add to list
     3139            spare[numberRemaining]=alpha;
     3140            index[numberRemaining++]=iSequence+addSequence;
     3141          }
     3142          break;
     3143        }
    30633144      }
    30643145    }
     
    31563237  double badFree=0.0;
    31573238  alpha_=0.0;
    3158   if (spareIntArray_[0]!=-1) {
     3239  if (spareIntArray_[0]>=0) {
    31593240    numberRemaining = dualColumn0(rowArray,columnArray,spareArray,
    31603241                                  acceptablePivot,upperTheta,bestPossible,badFree);
     
    31653246    upperTheta = spareDoubleArray_[0];
    31663247    bestPossible = spareDoubleArray_[1];
    3167     theta_ = spareDoubleArray_[2];
    3168     alpha_ = spareDoubleArray_[3];
    3169     sequenceIn_ = spareIntArray_[1];
     3248    if (spareIntArray_[0]==-1) {
     3249      theta_ = spareDoubleArray_[2];
     3250      alpha_ = spareDoubleArray_[3];
     3251      sequenceIn_ = spareIntArray_[1];
     3252    } else {
     3253#if 0
     3254      int n=numberRemaining;
     3255      double u=upperTheta;
     3256      double b=bestPossible;
     3257      upperTheta = 1.0e31;
     3258      bestPossible=0.0;
     3259      numberRemaining = dualColumn0(rowArray,columnArray,spareArray,
     3260                                    acceptablePivot,upperTheta,bestPossible,badFree);
     3261      assert (n==numberRemaining);
     3262      assert (fabs(b-bestPossible)<1.0e-7);
     3263      assert (fabs(u-upperTheta)<1.0e-7);
     3264#endif
     3265    }
    31703266  }
    31713267  // switch off
     
    40954191    factorization_->pivotTolerance(newTolerance);
    40964192  }
     4193  bestObjectiveValue_=CoinMax(bestObjectiveValue_,
     4194                              objectiveValue_-bestPossibleImprovement_);
    40974195  bool reallyBadProblems=false;
    40984196  // Double check infeasibility if no action
     
    43364434    handler_->message()<<CoinMessageEol;
    43374435  }
     4436#if 0
     4437  printf("IT %d %g %g(%d) %g(%d)\n",
     4438         numberIterations_,objectiveValue(),
     4439         sumPrimalInfeasibilities_,numberPrimalInfeasibilities_,
     4440         sumDualInfeasibilities_,numberDualInfeasibilities_);
     4441#endif   
    43384442#ifdef CLP_REPORT_PROGRESS
    43394443    if (ixxxxxx>=ixxyyyy-4&&ixxxxxx<=ixxyyyy) {
     
    58075911int ClpSimplexDual::fastDual(bool alwaysFinish)
    58085912{
     5913  bestObjectiveValue_ = objectiveValue_;
    58095914  algorithm_ = -1;
    58105915  secondaryStatus_=0;
     
    58275932  double saveDualBound = dualBound_;
    58285933
     5934  // Start can skip some things in transposeTimes
     5935  specialOptions_ |= 131072;
    58295936  if (alphaAccuracy_!=-1.0)
    58305937    alphaAccuracy_ = 1.0;
     
    60196126  dontFactorizePivots_ = saveDont;
    60206127  dualBound_ = saveDualBound;
     6128  // Stop can skip some things in transposeTimes
     6129  specialOptions_ &= ~131072;
     6130  if (problemStatus_==3)
     6131    objectiveValue_ = CoinMax(bestObjectiveValue_,objectiveValue_-bestPossibleImprovement_);
    60216132  return returnCode;
    60226133}
  • trunk/Clp/src/ClpSimplexOther.cpp

    r1370 r1376  
    13001300    small = new ClpSimplex(this,numberRows2,whichRow,
    13011301                     numberColumns2,whichColumn,true,false);
     1302#if 0
     1303    ClpPackedMatrix * rowCopy = dynamic_cast<ClpPackedMatrix *>(rowCopy_);
     1304    if (rowCopy) {
     1305      assert(!small->rowCopy());
     1306      small->setNewRowCopy(new ClpPackedMatrix(*rowCopy,numberRows2,whichRow,
     1307                                              numberColumns2,whichColumn));
     1308    }
     1309#endif
    13021310    // Set some stuff
    13031311    small->setDualBound(dualBound_);
  • trunk/Clp/src/ClpSimplexPrimal.cpp

    r1370 r1376  
    222222  }
    223223#endif
     224  // Start can skip some things in transposeTimes
     225  specialOptions_ |= 131072;
    224226  if (!startup(ifValuesPass,startFinishOptions)) {
    225227   
     
    483485    computeDuals(NULL);
    484486  }
     487  // Stop can skip some things in transposeTimes
     488  specialOptions_ &= ~131072;
    485489  // clean up
    486490  unflag();
     
    22402244    }
    22412245#ifdef KEEP_GOING_IF_FIXED
    2242     printf("ratio number diff rhs %g (%d %d fixed), element ratio %g\n",((double)number)/((double) numberRows_),
    2243            numberRows_,nFixed,elementRatio);
     2246    //printf("ratio number diff rhs %g (%d %d fixed), element ratio %g\n",((double)number)/((double) numberRows_),
     2247    //   numberRows_,nFixed,elementRatio);
    22442248#endif
    22452249    if (number*4>numberRows_||elementRatio>1.0e12) {
     
    22792283        last=sort[i];
    22802284      }
    2281       printf("cratio number diff bounds %g (%d %d fixed)\n",((double)number)/((double) numberColumns_),
    2282              numberColumns_,nFixed);
     2285      //printf("cratio number diff bounds %g (%d %d fixed)\n",((double)number)/((double) numberColumns_),
     2286      //     numberColumns_,nFixed);
    22832287    }
    22842288#endif
     
    32693273#if 0
    32703274  // if so - put in any superbasic costed slacks
     3275  // Start can skip some things in transposeTimes
     3276  specialOptions_ |= 131072;
    32713277  if (ifValuesPass&&specialOptions_<0x01000000) {
    32723278    // Get column copy
     
    36083614    computeDuals(NULL);
    36093615  }
     3616  // Stop can skip some things in transposeTimes
     3617  specialOptions_ &= ~131072;
    36103618  // clean up
    36113619  unflag();
Note: See TracChangeset for help on using the changeset viewer.