Changeset 210 for branches


Ignore:
Timestamp:
Sep 19, 2003 5:01:35 PM (16 years ago)
Author:
forrest
Message:

Trying

Location:
branches/pre
Files:
1 added
28 edited

Legend:

Unmodified
Added
Removed
  • branches/pre/ClpDualRowSteepest.cpp

    r208 r210  
    285285              //largestWeight = dubious;
    286286            }
     287          } else {
     288            // just to make sure we don't exit before got something
     289            numberWanted++;
    287290          }
    288291        }
  • branches/pre/ClpMatrixBase.cpp

    r205 r210  
    107107  return weights;
    108108}
     109// Append Columns
     110void
     111ClpMatrixBase::appendCols(int number, const CoinPackedVectorBase * const * columns)
     112{
     113  std::cerr<<"appendCols not supported - ClpMatrixBase"<<std::endl;
     114  abort();
     115}
     116// Append Rows
     117void
     118ClpMatrixBase::appendRows(int number, const CoinPackedVectorBase * const * rows)
     119{
     120  std::cerr<<"appendRows not supported - ClpMatrixBase"<<std::endl;
     121  abort();
     122}
  • branches/pre/ClpMessage.cpp

    r208 r210  
    2828  {CLP_DUAL_CHECKB,12,2,"New dual bound of %g"},
    2929  {CLP_DUAL_ORIGINAL,13,3,"Going back to original objective"},
    30   {CLP_SIMPLEX_PERTURB,14,1,"Perturbing problem by %g %% of %g - largest change %g (%% %g) - largest zero change %g"},
     30  {CLP_SIMPLEX_PERTURB,14,1,"Perturbing problem by %g %% of %g - largest nonzero change %g (%% %g) - largest zero change %g"},
    3131  {CLP_PRIMAL_ORIGINAL,15,2,"Going back to original tolerance"},
    3232  {CLP_PRIMAL_WEIGHT,16,2,"New infeasibility weight of %g"},
     
    6868  {CLP_TIMING,32,1,"%s objective %.9g - %d iterations time %.2f2%?, Presolve %.2f%?, Idiot %.2f%?"},
    6969  {CLP_INTERVAL_TIMING,33,2,"%s took %.2f seconds (total %.2f)"},
     70  {CLP_SPRINT,34,1,"Pass %d took %d iterations, objective %g, dual infeasibilities %g( %d)"},
    7071  {CLP_DUMMY_END,999999,0,""}
    7172};
  • branches/pre/ClpModel.cpp

    r208 r210  
    1313
    1414#include "CoinHelperFunctions.hpp"
     15#include "CoinTime.hpp"
    1516#include "ClpModel.hpp"
    1617#include "ClpPackedMatrix.hpp"
     
    4950  solveType_(0),
    5051  problemStatus_(-1),
     52  secondaryStatus_(0),
    5153  lengthNames_(0),
    5254  defaultHandler_(true),
     
    319321  solveType_ = rhs.solveType_;
    320322  problemStatus_ = rhs.problemStatus_;
     323  secondaryStatus_ = rhs.secondaryStatus_;
    321324  numberRows_ = rhs.numberRows_;
    322325  numberColumns_ = rhs.numberColumns_;
     
    361364    status_ = ClpCopyOfArray( rhs.status_,numberColumns_+numberRows_);
    362365    ray_ = NULL;
    363     if (problemStatus_==1)
     366    if (problemStatus_==1&&!secondaryStatus_)
    364367      ray_ = ClpCopyOfArray (rhs.ray_,numberRows_);
    365368    else if (problemStatus_==2)
     
    420423  otherModel.numberIterations_ = numberIterations_;
    421424  otherModel.problemStatus_ = problemStatus_;
     425  otherModel.secondaryStatus_ = secondaryStatus_;
    422426  rowActivity_ = NULL;
    423427  columnActivity_ = NULL;
     
    648652    // set state back to unknown
    649653    problemStatus_ = -1;
     654    secondaryStatus_ = 0;
    650655    delete [] ray_;
    651656    ray_ = NULL;
     
    715720  // set state back to unknown
    716721  problemStatus_ = -1;
     722  secondaryStatus_ = 0;
    717723  delete [] ray_;
    718724  ray_ = NULL;
     
    756762  // set state back to unknown
    757763  problemStatus_ = -1;
     764  secondaryStatus_ = 0;
    758765  delete [] ray_;
    759766  ray_ = NULL;
     
    856863  if (!matrix_)
    857864    createEmptyMatrix();
    858   // Use matrix() to get round virtual problem
    859   matrix()->appendRows(number,rows);
     865  matrix_->appendRows(number,rows);
    860866}
    861867// Add columns
     
    967973  if (!matrix_)
    968974    createEmptyMatrix();
    969   // Use matrix() to get round virtual problem
    970   matrix()->appendCols(number,columns);
     975  matrix_->appendCols(number,columns);
    971976}
    972977// Infeasibility/unbounded ray (NULL returned if none/wrong)
     
    975980{
    976981  double * array = NULL;
    977   if (problemStatus_==1)
     982  if (problemStatus_==1&&!secondaryStatus_)
    978983    array = ClpCopyOfArray(ray_,numberRows_);
    979984  return array;
     
    10481053  }
    10491054  CoinMpsIO m;
     1055  bool savePrefix =m.messageHandler()->prefix();
     1056  m.messageHandler()->setPrefix(handler_->prefix());
    10501057  double time1 = CoinCpuTime(),time2;
    10511058  int status=m.readMps(fileName,"");
     1059  m.messageHandler()->setPrefix(savePrefix);
    10521060  if (!status||ignoreErrors) {
    10531061    loadProblem(*m.getMatrixByCol(),
     
    13341342  solveType_ = rhs->solveType_;
    13351343  problemStatus_ = rhs->problemStatus_;
     1344  secondaryStatus_ = rhs->secondaryStatus_;
    13361345  // check valid lists
    13371346  int numberBad=0;
     
    14081417  delete [] columnStatus;
    14091418  ray_ = NULL;
    1410   if (problemStatus_==1)
     1419  if (problemStatus_==1&&!secondaryStatus_)
    14111420    ray_ = whichDouble (rhs->ray_,numberRows,whichRow);
    14121421  else if (problemStatus_==2)
  • branches/pre/ClpNonLinearCost.cpp

    r192 r210  
    440440        }
    441441      } else {
    442         if (fabs(value-upperValue)<=fabs(value-lowerValue)) {
    443           solution[iSequence] = upperValue;
    444         } else {
    445           model_->setStatus(iSequence,ClpSimplex::atLowerBound);
    446           solution[iSequence] = lowerValue;
    447         }
     442        // Set to nearest and make at upper bound
     443        int kRange;
     444        iRange=-1;
     445        double nearest = COIN_DBL_MAX;
     446        for (kRange=start; kRange<end;kRange++) {
     447          if (fabs(lower_[kRange]-value)<nearest) {
     448            nearest = fabs(lower_[kRange]-value);
     449            iRange=kRange;
     450          }
     451        }
     452        assert (iRange>=0);
     453        iRange--;
     454        solution[iSequence]=lower_[iRange+1];
    448455      }
    449456      break;
     
    460467        }
    461468      } else {
    462         if (fabs(value-lowerValue)<=fabs(value-upperValue)) {
    463           solution[iSequence] = lowerValue;
    464         } else {
    465           model_->setStatus(iSequence,ClpSimplex::atUpperBound);
    466           solution[iSequence] = upperValue;
    467         }
     469        // Set to nearest and make at lower bound
     470        int kRange;
     471        iRange=-1;
     472        double nearest = COIN_DBL_MAX;
     473        for (kRange=start; kRange<end;kRange++) {
     474          if (fabs(lower_[kRange]-value)<nearest) {
     475            nearest = fabs(lower_[kRange]-value);
     476            iRange=kRange;
     477          }
     478        }
     479        assert (iRange>=0);
     480        solution[iSequence]=lower_[iRange];
    468481      }
    469482      break;
    470483    case ClpSimplex::isFixed:
     484      if (toNearest) {
     485        // Set to true fixed
     486        for (iRange=start; iRange<end;iRange++) {
     487          if (lower_[iRange]==lower_[iRange+1])
     488            break;
     489        }
     490        solution[iSequence]=lower_[iRange];
     491      }
    471492      break;
    472493    }
     
    745766  return difference;
    746767}
     768/* Sets bounds and cost for outgoing variable
     769   may change value
     770   Returns direction */
     771int
     772ClpNonLinearCost::setOneOutgoing(int iPivot, double & value)
     773{
     774  assert (model_!=NULL);
     775  double primalTolerance = model_->currentPrimalTolerance();
     776  // get where in bound sequence
     777  int iRange;
     778  int currentRange = whichRange_[iPivot];
     779  int start = start_[iPivot];
     780  int end = start_[iPivot+1]-1;
     781  // Set perceived direction out
     782  int direction;
     783  if (value<=lower_[currentRange]+1.001*primalTolerance) {
     784    direction=1;
     785  } else if (value>=lower_[currentRange+1]-1.001*primalTolerance) {
     786    direction=-1;
     787  } else {
     788    // odd
     789    direction=0;
     790  }
     791  // If fixed try and get feasible
     792  if (lower_[start+1]==lower_[start+2]&&fabs(value-lower_[start+1])<1.001*primalTolerance) {
     793    iRange =start+1;
     794  } else {
     795    // See if exact
     796    for (iRange=start; iRange<end;iRange++) {
     797      if (value==lower_[iRange+1]) {
     798        // put in better range
     799        if (infeasible(iRange)&&iRange==start)
     800          iRange++;
     801        break;
     802      }
     803    }
     804    if (iRange==end) {
     805      // not exact
     806      for (iRange=start; iRange<end;iRange++) {
     807        if (value<=lower_[iRange+1]+primalTolerance) {
     808          // put in better range
     809          if (value>=lower_[iRange+1]-primalTolerance&&infeasible(iRange)&&iRange==start)
     810            iRange++;
     811          break;
     812        }
     813      }
     814    }
     815  }
     816  assert(iRange<end);
     817  whichRange_[iPivot]=iRange;
     818  if (iRange!=currentRange) {
     819    if (infeasible(iRange))
     820      numberInfeasibilities_++;
     821    if (infeasible(currentRange))
     822      numberInfeasibilities_--;
     823  }
     824  double & lower = model_->lowerAddress(iPivot);
     825  double & upper = model_->upperAddress(iPivot);
     826  double & cost = model_->costAddress(iPivot);
     827  lower = lower_[iRange];
     828  upper = lower_[iRange+1];
     829  if (upper==lower) {
     830    value=upper;
     831  } else {
     832    // set correctly
     833    if (fabs(value-lower)<=primalTolerance*1.001){
     834      value = min(value,lower+primalTolerance);
     835    } else if (fabs(value-upper)<=primalTolerance*1.001){
     836      value = max(value,upper-primalTolerance);
     837    } else {
     838      printf("*** variable wandered off bound %g %g %g!\n",
     839             lower,value,upper);
     840      if (value-lower<=upper-value)
     841        value = lower+primalTolerance;
     842      else
     843        value = upper-primalTolerance;
     844    }
     845  }
     846  double difference = cost-cost_[iRange];
     847  cost = cost_[iRange];
     848  changeCost_ += value*difference;
     849  return direction;
     850}
    747851// Returns nearest bound
    748852double
  • branches/pre/ClpPlusMinusOneMatrix.cpp

    r208 r210  
    119119    delete [] indices_;
    120120    // put in message
    121     printf("Not all +-1 - can test if indices_ null\n");
    122121    indices_=NULL;
    123122    numberRows_=0;
     
    275274        // bad column
    276275        numberBad++;
     276        printf("%d %d %d %d\n",iColumn,numberMajor,numberMajor1,kColumn);
    277277      }
    278278    }
     
    471471    else if (numberRows*2<numberColumns)
    472472      factor=0.2;
    473     if (model->numberIterations()%50==0)
    474       printf("%d nonzero\n",numberInRowArray);
    475473  }
    476474  if (numberInRowArray>factor*numberRows||!rowCopy) {
     
    11191117ClpPlusMinusOneMatrix::deleteCols(const int numDel, const int * indDel)
    11201118{
    1121   abort();
     1119  int iColumn;
     1120  int newSize=startPositive_[numberColumns_];;
     1121  int numberBad=0;
     1122  // Use array to make sure we can have duplicates
     1123  int * which = new int[numberColumns_];
     1124  memset(which,0,numberColumns_*sizeof(int));
     1125  int nDuplicate=0;
     1126  for (iColumn=0;iColumn<numDel;iColumn++) {
     1127    int jColumn = indDel[iColumn];
     1128    if (jColumn<0||jColumn>=numberColumns_) {
     1129      numberBad++;
     1130    } else {
     1131      newSize -= startPositive_[jColumn+1]-startPositive_[jColumn];
     1132      if (which[jColumn])
     1133        nDuplicate++;
     1134      else
     1135        which[jColumn]=1;
     1136    }
     1137  }
     1138  if (numberBad)
     1139    throw CoinError("Indices out of range", "deleteCols", "ClpPlusMinusOneMatrix");
     1140  int newNumber = numberColumns_-numDel+nDuplicate;
     1141  // Get rid of temporary arrays
     1142  delete [] lengths_;
     1143  lengths_=NULL;
     1144  delete [] elements_;
     1145  elements_= NULL;
     1146  int * newPositive = new int [newNumber+1];
     1147  int * newNegative = new int [newNumber];
     1148  int * newIndices = new int [newSize];
     1149  newNumber=0;
     1150  newSize=0;
     1151  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
     1152    if (!which[iColumn]) {
     1153      int start,end;
     1154      int i;
     1155      start = startPositive_[iColumn];
     1156      end=startNegative_[iColumn];
     1157      newPositive[newNumber]=newSize;
     1158      for (i=start;i<end;i++)
     1159        newIndices[newSize++]=indices_[i];
     1160      start = startNegative_[iColumn];
     1161      end=startPositive_[iColumn+1];
     1162      newNegative[newNumber++]=newSize;
     1163      for (i=start;i<end;i++)
     1164        newIndices[newSize++]=indices_[i];
     1165    }
     1166  }
     1167  newPositive[newNumber]=newSize;
     1168  delete [] which;
     1169  delete [] startPositive_;
     1170  startPositive_= newPositive;
     1171  delete [] startNegative_;
     1172  startNegative_= newNegative;
     1173  delete [] indices_;
     1174  indices_= newIndices;
     1175  numberColumns_ = newNumber;
    11221176}
    11231177/* Delete the rows whose indices are listed in <code>indDel</code>. */
     
    11251179ClpPlusMinusOneMatrix::deleteRows(const int numDel, const int * indDel)
    11261180{
    1127   abort();
     1181  int iRow;
     1182  int numberBad=0;
     1183  // Use array to make sure we can have duplicates
     1184  int * which = new int[numberRows_];
     1185  memset(which,0,numberRows_*sizeof(int));
     1186  int nDuplicate=0;
     1187  for (iRow=0;iRow<numDel;iRow++) {
     1188    int jRow = indDel[iRow];
     1189    if (jRow<0||jRow>=numberRows_) {
     1190      numberBad++;
     1191    } else {
     1192      if (which[jRow])
     1193        nDuplicate++;
     1194      else
     1195        which[jRow]=1;
     1196    }
     1197  }
     1198  if (numberBad)
     1199    throw CoinError("Indices out of range", "deleteCols", "ClpPlusMinusOneMatrix");
     1200  int iElement;
     1201  int numberElements=startPositive_[numberColumns_];
     1202  int newSize=0;
     1203  for (iElement=0;iElement<numberElements;iElement++) {
     1204    iRow = indices_[iElement];
     1205    if (!which[iRow])
     1206      newSize++;
     1207  }
     1208  int newNumber = numberRows_-numDel+nDuplicate;
     1209  // Get rid of temporary arrays
     1210  delete [] lengths_;
     1211  lengths_=NULL;
     1212  delete [] elements_;
     1213  elements_= NULL;
     1214  int * newIndices = new int [newSize];
     1215  newSize=0;
     1216  int iColumn;
     1217  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
     1218    int start,end;
     1219    int i;
     1220    start = startPositive_[iColumn];
     1221    end=startNegative_[iColumn];
     1222    startPositive_[newNumber]=newSize;
     1223    for (i=start;i<end;i++) {
     1224        iRow = indices_[i];
     1225        if (!which[iRow])
     1226          newIndices[newSize++]=iRow;
     1227    }
     1228    start = startNegative_[iColumn];
     1229    end=startPositive_[iColumn+1];
     1230    startNegative_[newNumber]=newSize;
     1231    for (i=start;i<end;i++) {
     1232        iRow = indices_[i];
     1233        if (!which[iRow])
     1234          newIndices[newSize++]=iRow;
     1235    }
     1236  }
     1237  startPositive_[numberColumns_]=newSize;
     1238  delete [] which;
     1239  delete [] indices_;
     1240  indices_= newIndices;
     1241  numberRows_ = newNumber;
    11281242}
    11291243bool
     
    11811295  return weights;
    11821296}
     1297// Append Columns
     1298void
     1299ClpPlusMinusOneMatrix::appendCols(int number, const CoinPackedVectorBase * const * columns)
     1300{
     1301  int iColumn;
     1302  int size=0;
     1303  int numberBad=0;
     1304  for (iColumn=0;iColumn<number;iColumn++) {
     1305    int n=columns[iColumn]->getNumElements();
     1306    const double * element = columns[iColumn]->getElements();
     1307    size += n;
     1308    int i;
     1309    for (i=0;i<n;i++) {
     1310      if (fabs(element[i])!=1.0)
     1311        numberBad++;
     1312    }
     1313  }
     1314  if (numberBad)
     1315    throw CoinError("Not +- 1", "appendCols", "ClpPlusMinusOneMatrix");
     1316  // Get rid of temporary arrays
     1317  delete [] lengths_;
     1318  lengths_=NULL;
     1319  delete [] elements_;
     1320  elements_= NULL;
     1321  int numberNow = startPositive_[numberColumns_];
     1322  int * temp;
     1323  temp = new int [numberColumns_+1+number];
     1324  memcpy(temp,startPositive_,(numberColumns_+1)*sizeof(int));
     1325  delete [] startPositive_;
     1326  startPositive_= temp;
     1327  temp = new int [numberColumns_+number];
     1328  memcpy(temp,startNegative_,numberColumns_*sizeof(int));
     1329  delete [] startNegative_;
     1330  startNegative_= temp;
     1331  temp = new int [numberNow+size];
     1332  memcpy(temp,indices_,numberNow*sizeof(int));
     1333  delete [] indices_;
     1334  indices_= temp;
     1335  // now add
     1336  size=numberNow;
     1337  for (iColumn=0;iColumn<number;iColumn++) {
     1338    int n=columns[iColumn]->getNumElements();
     1339    const int * row = columns[iColumn]->getIndices();
     1340    const double * element = columns[iColumn]->getElements();
     1341    int i;
     1342    for (i=0;i<n;i++) {
     1343      if (element[i]==1.0)
     1344        indices_[size++]=row[i];
     1345    }
     1346    startNegative_[iColumn+numberColumns_]=size;
     1347    for (i=0;i<n;i++) {
     1348      if (element[i]==-1.0)
     1349        indices_[size++]=row[i];
     1350    }
     1351    startPositive_[iColumn+numberColumns_+1]=size;
     1352  }
     1353 
     1354  numberColumns_ += number;
     1355}
     1356// Append Rows
     1357void
     1358ClpPlusMinusOneMatrix::appendRows(int number, const CoinPackedVectorBase * const * rows)
     1359{
     1360  // Allocate arrays to use for counting
     1361  int * countPositive = new int [numberColumns_+1];
     1362  memset(countPositive,0,numberColumns_*sizeof(int));
     1363  int * countNegative = new int [numberColumns_];
     1364  memset(countNegative,0,numberColumns_*sizeof(int));
     1365  int iRow;
     1366  int size=0;
     1367  int numberBad=0;
     1368  for (iRow=0;iRow<number;iRow++) {
     1369    int n=rows[iRow]->getNumElements();
     1370    const int * column = rows[iRow]->getIndices();
     1371    const double * element = rows[iRow]->getElements();
     1372    size += n;
     1373    int i;
     1374    for (i=0;i<n;i++) {
     1375      int iColumn = column[i];
     1376      if (element[i]==1.0)
     1377        countPositive[iColumn++];
     1378      else if (element[i]==-1.0)
     1379        countNegative[iColumn++];
     1380      else
     1381        numberBad++;
     1382    }
     1383  }
     1384  if (numberBad)
     1385    throw CoinError("Not +- 1", "appendRows", "ClpPlusMinusOneMatrix");
     1386  // Get rid of temporary arrays
     1387  delete [] lengths_;
     1388  lengths_=NULL;
     1389  delete [] elements_;
     1390  elements_= NULL;
     1391  int numberNow = startPositive_[numberColumns_];
     1392  int * newIndices = new int [numberNow+size];
     1393  // Update starts and turn counts into positions
     1394  // also move current indices
     1395  int iColumn;
     1396  int numberAdded=0;
     1397  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
     1398    int n,now,move;
     1399    now = startPositive_[iColumn];
     1400    move = startNegative_[iColumn]-now;
     1401    n = countPositive[iColumn];
     1402    startPositive_[iColumn] += numberAdded;
     1403    memcpy(indices_+now,newIndices+startPositive_[iColumn],move*sizeof(int));
     1404    countPositive[iColumn]= startNegative_[iColumn]+numberAdded;
     1405    numberAdded += n;
     1406    now = startNegative_[iColumn];
     1407    move = startPositive_[iColumn+1]-now;
     1408    n = countNegative[iColumn];
     1409    startNegative_[iColumn] += numberAdded;
     1410    memcpy(indices_+now,newIndices+startNegative_[iColumn],move*sizeof(int));
     1411    countNegative[iColumn]= startPositive_[iColumn+1]+numberAdded;
     1412    numberAdded += n;
     1413  }
     1414  delete [] indices_;
     1415  indices_ = newIndices;
     1416  startPositive_[numberColumns_] += numberAdded;
     1417  // Now put in
     1418  for (iRow=0;iRow<number;iRow++) {
     1419    int newRow = numberRows_+iRow;
     1420    int n=rows[iRow]->getNumElements();
     1421    const int * column = rows[iRow]->getIndices();
     1422    const double * element = rows[iRow]->getElements();
     1423    int i;
     1424    for (i=0;i<n;i++) {
     1425      int iColumn = column[i];
     1426      int put;
     1427      if (element[i]==1.0) {
     1428        put = countPositive[iColumn];
     1429        countPositive[iColumn] = put+1;
     1430      } else {
     1431        put = countNegative[iColumn];
     1432        countNegative[iColumn] = put+1;
     1433      }
     1434      indices_[put]=newRow;
     1435    }
     1436  }
     1437  delete [] countPositive;
     1438  delete [] countNegative;
     1439  numberRows_ += number;
     1440}
  • branches/pre/ClpPrimalColumnSteepest.cpp

    r197 r210  
    3535    state_(-1),
    3636    mode_(mode),
     37    numberSwitched_(0),
    3738    pivotSequence_(-1),
    3839    savedPivotSequence_(-1),
     
    5051  state_=rhs.state_;
    5152  mode_ = rhs.mode_;
     53  numberSwitched_ = rhs.numberSwitched_;
    5254  model_ = rhs.model_;
    5355  pivotSequence_ = rhs.pivotSequence_;
     
    6870    savedWeights_= new double[number];
    6971    ClpDisjointCopyN(rhs.savedWeights_,number,savedWeights_);
    70     if (!mode_) {
     72    if (mode_!=1) {
    7173      reference_ = new unsigned int[(number+31)>>5];
    7274      memcpy(reference_,rhs.reference_,((number+31)>>5)*sizeof(unsigned int));
     
    105107    state_=rhs.state_;
    106108    mode_ = rhs.mode_;
     109    numberSwitched_ = rhs.numberSwitched_;
    107110    model_ = rhs.model_;
    108111    pivotSequence_ = rhs.pivotSequence_;
     
    129132      savedWeights_= new double[number];
    130133      ClpDisjointCopyN(rhs.savedWeights_,number,savedWeights_);
    131       if (!mode_) {
     134      if (mode_!=1) {
    132135        reference_ = new unsigned int[(number+31)>>5];
    133136        memcpy(reference_,rhs.reference_,
     
    164167  // dj could be very small (or even zero - take care)
    165168  double dj = model_->dualIn();
    166   double * infeas = infeasible_->denseVector();
    167169  double tolerance=model_->currentDualTolerance();
    168170  // we can't really trust infeasibilities if there is dual error
     
    172174  tolerance = tolerance +  error;
    173175  int pivotRow = model_->pivotRow();
    174 
     176  // Local copy of mode so can decide what to do
     177  int switchType = mode_;
     178  if (mode_==4) {
     179    if (numberSwitched_) {
     180      switchType=3;
     181    } else {
     182      // See if to switch
     183      int numberElements = model_->factorization()->numberElements();
     184      int numberRows = model_->numberRows();
     185      // ratio is done on number of columns here
     186      int numberColumns = model_->numberColumns();
     187      double ratio = (double) numberElements/(double) numberColumns;
     188      //double ratio = (double) numberElements/(double) model_->clpMatrix()->getNumElements();
     189      int numberWanted;
     190      bool doPartial=true;
     191      int numberTotal = numberColumns+model_->numberRows();
     192      if (ratio<0.1) {
     193        numberWanted = max(100,numberTotal/200);
     194      } else if (ratio<1.0) {
     195        numberWanted = max(100,numberTotal/20);
     196      } else if (ratio<2.0) {
     197        numberWanted = max(200,numberTotal/10);
     198      } else {
     199        doPartial=false; // switch off
     200        numberWanted=0; // just for compliler message
     201      }
     202      if (doPartial) {
     203        if (model_->numberIterations()%100==0)
     204          printf("numels %d ratio %g wanted %d\n",numberElements,ratio,numberWanted);
     205        // No do partial
     206        // Always do slacks
     207        int anyUpdates;
     208        double * dual = model_->dualRowSolution();
     209        if (updates->getNumElements()) {
     210          // would have to have two goes for devex, three for steepest
     211          anyUpdates=2;
     212          // add in pivot contribution
     213          if (pivotRow>=0)
     214            updates->add(pivotRow,-dj);
     215        } else if (pivotRow>=0) {
     216          if (fabs(dj)>1.0e-15) {
     217            // some dj
     218            updates->insert(pivotRow,-dj);
     219            if (fabs(dj)>1.0e-6) {
     220              // reasonable size
     221              anyUpdates=1;
     222            } else {
     223              // too small
     224              anyUpdates=2;
     225            }
     226          } else {
     227            // zero dj
     228            anyUpdates=-1;
     229          }
     230        } else if (pivotSequence_>=0){
     231          // just after re-factorization
     232          anyUpdates=-1;
     233        } else {
     234          // sub flip - nothing to do
     235          anyUpdates=0;
     236        }
     237
     238        // For now (as we can always switch off)
     239        assert (!model_->nonLinearCost()->lookBothWays());
     240        double * slackDj =model_->djRegion(0);
     241        double * dj =model_->djRegion(1);
     242        if (anyUpdates>0) {
     243          model_->factorization()->updateColumnTranspose(spareRow2,updates);
     244   
     245          number = updates->getNumElements();
     246          index = updates->getIndices();
     247          updateBy = updates->denseVector();
     248         
     249
     250          for (j=0;j<number;j++) {
     251            int iSequence = index[j];
     252            double value = updateBy[iSequence];
     253            slackDj[iSequence] -= value;
     254            dual[iSequence] -= value;
     255            updateBy[iSequence]=0.0;
     256          }
     257          updates->setNumElements(0);
     258        }
     259       
     260        int iPass;
     261        // Setup two passes AND treat slacks differently
     262        int start[2][4];
     263        // slacks
     264        start[0][1]=numberRows;
     265        start[0][2]=0;
     266        double dstart = ((double) numberRows) * CoinDrand48();
     267        start[0][0]=(int) dstart;
     268        start[0][3]=start[0][0];
     269        for (iPass=0;iPass<4;iPass++)
     270          start[0][iPass] += numberColumns;
     271        start[1][1]=numberColumns;
     272        start[1][2]=0;
     273        dstart = ((double) numberColumns) * CoinDrand48();
     274        start[1][0]=(int) dstart;
     275        start[1][3]=start[1][0];
     276        int bestSequence=-1;
     277        double bestDj=tolerance;
     278        // We are going to fake it so pi looks sparse
     279        double * saveDense = updates->denseVector();
     280        updates->setDenseVector(dual);
     281        int kChunk = max(2000,numberWanted);
     282        for (iPass=0;iPass<2;iPass++) {
     283          // Slacks
     284          int end = start[0][2*iPass+1];
     285          for (int iSequence=start[0][2*iPass];iSequence<end;iSequence++) {
     286            double value;
     287            ClpSimplex::Status status = model_->getStatus(iSequence);
     288         
     289            switch(status) {
     290             
     291            case ClpSimplex::basic:
     292            case ClpSimplex::isFixed:
     293              break;
     294            case ClpSimplex::isFree:
     295            case ClpSimplex::superBasic:
     296              value = dj[iSequence];
     297              if (fabs(value)>FREE_ACCEPT*tolerance) {
     298                numberWanted--;
     299                // we are going to bias towards free (but only if reasonable)
     300                value = fabs(value*FREE_BIAS);
     301                if (value>bestDj) {
     302                  // check flagged variable and correct dj
     303                  if (!model_->flagged(iSequence)) {
     304                    bestDj=value;
     305                    bestSequence = iSequence;
     306                  } else {
     307                    // just to make sure we don't exit before got something
     308                    numberWanted++;
     309                  }
     310                }
     311              }
     312              break;
     313            case ClpSimplex::atUpperBound:
     314              value = dj[iSequence];
     315              if (value>tolerance) {
     316                numberWanted--;
     317                if (value>bestDj) {
     318                  // check flagged variable and correct dj
     319                  if (!model_->flagged(iSequence)) {
     320                    bestDj=value;
     321                    bestSequence = iSequence;
     322                  } else {
     323                    // just to make sure we don't exit before got something
     324                    numberWanted++;
     325                  }
     326                }
     327              }
     328              break;
     329            case ClpSimplex::atLowerBound:
     330              value = dj[iSequence];
     331              if (value<-tolerance) {
     332                numberWanted--;
     333                if (-value>bestDj) {
     334                  // check flagged variable and correct dj
     335                  if (!model_->flagged(iSequence)) {
     336                    bestDj=-value;
     337                    bestSequence = iSequence;
     338                  } else {
     339                    // just to make sure we don't exit before got something
     340                    numberWanted++;
     341                  }
     342                }
     343              }
     344            }
     345            if (!numberWanted)
     346              break;
     347          }
     348          if (!numberWanted)
     349            break;
     350          // For Columns break up into manageable chunks
     351          end = start[1][2*iPass+1];
     352          index = spareColumn1->getIndices();
     353          updateBy = spareColumn1->denseVector();
     354          for (int jSequence=start[1][2*iPass];jSequence<end;jSequence+=kChunk) {
     355            int endThis = min(end,jSequence+kChunk);
     356            number=0;
     357            for (int iSequence=jSequence;iSequence<endThis;iSequence++) {
     358              ClpSimplex::Status status = model_->getStatus(iSequence);
     359         
     360              if(status!=ClpSimplex::basic&&status!=ClpSimplex::isFixed) {
     361                index[number++]=iSequence;
     362              }
     363            }
     364            spareColumn1->setNumElements(number);
     365            // get subset which have nonzero djs
     366            // updates is actually slack djs
     367            model_->clpMatrix()->subsetTransposeTimes(model_,updates,
     368                                                      spareColumn1,
     369                                                      spareColumn2);
     370            double * updateBy2 = spareColumn2->denseVector();
     371            const double * cost = model_->costRegion();
     372            for (int j=0;j<number;j++) {
     373              int iSequence = index[j];
     374              double value = cost[iSequence]-updateBy2[iSequence];
     375              updateBy2[iSequence]=0.0;
     376              ClpSimplex::Status status = model_->getStatus(iSequence);
     377         
     378              switch(status) {
     379               
     380              case ClpSimplex::basic:
     381              case ClpSimplex::isFixed:
     382                break;
     383              case ClpSimplex::isFree:
     384              case ClpSimplex::superBasic:
     385                if (fabs(value)>FREE_ACCEPT*tolerance) {
     386                  numberWanted--;
     387                  // we are going to bias towards free (but only if reasonable)
     388                  if (fabs(value*FREE_BIAS)>bestDj) {
     389                    // check flagged variable and correct dj
     390                    if (!model_->flagged(iSequence)) {
     391                      bestDj=value*FREE_BIAS;
     392                      dj[iSequence]=value;
     393                      bestSequence = iSequence;
     394                    } else {
     395                      // just to make sure we don't exit before got something
     396                      numberWanted++;
     397                    }
     398                  }
     399                }
     400                break;
     401              case ClpSimplex::atUpperBound:
     402                if (value>tolerance) {
     403                  numberWanted--;
     404                  if (value>bestDj) {
     405                    // check flagged variable and correct dj
     406                    if (!model_->flagged(iSequence)) {
     407                      bestDj=value;
     408                      dj[iSequence]=value;
     409                      bestSequence = iSequence;
     410                    } else {
     411                      // just to make sure we don't exit before got something
     412                      numberWanted++;
     413                    }
     414                  }
     415                }
     416                break;
     417              case ClpSimplex::atLowerBound:
     418                if (value<-tolerance) {
     419                  numberWanted--;
     420                  if (-value>bestDj) {
     421                    // check flagged variable and correct dj
     422                    if (!model_->flagged(iSequence)) {
     423                      bestDj=-value;
     424                      dj[iSequence]=value;
     425                      bestSequence = iSequence;
     426                    } else {
     427                      // just to make sure we don't exit before got something
     428                      numberWanted++;
     429                    }
     430                  }
     431                }
     432                break;
     433              }
     434              if (!numberWanted) {
     435                // erase rest
     436                for (;j<number;j++) {
     437                  int iSequence = index[j];
     438                  updateBy2[iSequence]=0.0;
     439                }
     440                break;
     441              }
     442            }
     443            if (!numberWanted)
     444              break;
     445          }
     446          if (!numberWanted)
     447            break;
     448        }
     449        spareColumn1->setNumElements(0);
     450        spareColumn2->setNumElements(0);
     451        updates->setDenseVector(saveDense);
     452        /*if (model_->numberIterations()%100==0)
     453          printf("%d best %g\n",bestSequence,bestDj);*/
     454       
     455#ifdef CLP_DEBUG
     456        if (bestSequence>=0) {
     457          if (model_->getStatus(bestSequence)==ClpSimplex::atLowerBound)
     458            assert(model_->reducedCost(bestSequence)<0.0);
     459          if (model_->getStatus(bestSequence)==ClpSimplex::atUpperBound)
     460            assert(model_->reducedCost(bestSequence)>0.0);
     461        }
     462#endif
     463        return bestSequence;
     464      } else {
     465        // Yes initialize
     466        pivotSequence_=-1;
     467        pivotRow=-1;
     468        numberSwitched_++;
     469        // Redo dualsolution
     470        model_->computeDuals(NULL);
     471        saveWeights(model_,4);
     472        printf("switching %d nel ratio %g\n",numberElements,ratio);
     473        updates->clear();
     474      }
     475    }
     476  }
    175477  int anyUpdates;
     478  double * infeas = infeasible_->denseVector();
    176479
    177480  if (updates->getNumElements()) {
     
    562865    // and we can see if reference
    563866    double referenceIn=0.0;
    564     if (!mode_&&reference(sequenceIn))
     867    if (switchType!=1&&reference(sequenceIn))
    565868      referenceIn=1.0;
    566869    // save outgoing weight round update
     
    604907      thisWeight += pivotSquared * devex_ + pivot * modification;
    605908      if (thisWeight<TRY_NORM) {
    606         if (mode_) {
     909        if (switchType==1) {
    607910          // steepest
    608911          thisWeight = max(TRY_NORM,ADD_ONE+pivotSquared);
     
    642945      thisWeight += pivotSquared * devex_ + pivot * modification;
    643946      if (thisWeight<TRY_NORM) {
    644         if (mode_) {
     947        if (switchType==1) {
    645948          // steepest
    646949          thisWeight = max(TRY_NORM,ADD_ONE+pivotSquared);
     
    7041007  }
    7051008  tolerance *= tolerance; // as we are using squares
    706   for (i=0;i<number;i++) {
    707     iSequence = index[i];
    708     double value = infeas[iSequence];
    709     double weight = weights_[iSequence];
    710     //weight=1.0;
    711     if (value>bestDj*weight&&value>tolerance) {
    712       // check flagged variable and correct dj
    713       if (!model_->flagged(iSequence)) {
    714         /*if (model_->numberIterations()%100==0)
    715           printf("chosen %g => %g\n",bestDj,value/weight);*/
    716         bestDj=value/weight;
    717         bestSequence = iSequence;
    718       }
    719     }
     1009
     1010  int numberWanted;
     1011  if (switchType<2 ) {
     1012    numberWanted = number+1;
     1013  } else if (switchType==2) {
     1014    numberWanted = max(2000,number/8);
     1015  } else {
     1016    int numberElements = model_->factorization()->numberElements();
     1017    double ratio = (double) numberElements/(double) model_->numberRows();
     1018    numberWanted = max(2000,number/8);
     1019    if (ratio<1.0) {
     1020      numberWanted = max(2000,number/20);
     1021    } else if (ratio>10.0) {
     1022      ratio = number * (ratio/80.0);
     1023      if (ratio>number)
     1024        numberWanted=number+1;
     1025      else
     1026        numberWanted = max(2000,(int) ratio);
     1027    }
     1028  }
     1029  int iPass;
     1030  // Setup two passes
     1031  int start[4];
     1032  start[1]=number;
     1033  start[2]=0;
     1034  double dstart = ((double) number) * CoinDrand48();
     1035  start[0]=(int) dstart;
     1036  start[3]=start[0];
     1037  //double largestWeight=0.0;
     1038  //double smallestWeight=1.0e100;
     1039  for (iPass=0;iPass<2;iPass++) {
     1040    int end = start[2*iPass+1];
     1041    for (i=start[2*iPass];i<end;i++) {
     1042      iSequence = index[i];
     1043      double value = infeas[iSequence];
     1044      if (value>tolerance) {
     1045        double weight = weights_[iSequence];
     1046        //weight=1.0;
     1047        if (value>bestDj*weight) {
     1048          // check flagged variable and correct dj
     1049          if (!model_->flagged(iSequence)) {
     1050            bestDj=value/weight;
     1051            bestSequence = iSequence;
     1052          } else {
     1053            // just to make sure we don't exit before got something
     1054            numberWanted++;
     1055          }
     1056        }
     1057      }
     1058      numberWanted--;
     1059      if (!numberWanted)
     1060        break;
     1061    }
     1062    if (!numberWanted)
     1063      break;
    7201064  }
    7211065  if (sequenceOut>=0) {
     
    7441088ClpPrimalColumnSteepest::saveWeights(ClpSimplex * model,int mode)
    7451089{
     1090  model_ = model;
     1091  if (mode_==4) {
     1092    if (mode==1&&!weights_)
     1093      numberSwitched_=0; // Reset
     1094    if(!numberSwitched_)
     1095      return;
     1096  }
    7461097  // alternateWeights_ is defined as indexed but is treated oddly
    7471098  // at times
    748   model_ = model;
    7491099  int numberRows = model_->numberRows();
    7501100  int numberColumns = model_->numberColumns();
     
    9281278ClpPrimalColumnSteepest::unrollWeights()
    9291279{
     1280  if (mode_==4&&!numberSwitched_)
     1281    return;
    9301282  double * saved = alternateWeights_->denseVector();
    9311283  int number = alternateWeights_->getNumElements();
     
    9541306ClpPrimalColumnSteepest::updateWeights(CoinIndexedVector * input)
    9551307{
     1308  // Local copy of mode so can decide what to do
     1309  int switchType = mode_;
     1310  if (mode_==4&&numberSwitched_)
     1311    switchType=3;
     1312  else if (mode_==4)
     1313    return;
    9561314  int number=input->getNumElements();
    9571315  int * which = input->getIndices();
     
    9711329
    9721330  if (pivotRow>=0) {
    973     if (mode_) {
     1331    if (switchType==1) {
    9741332      for (i=0;i<number;i++) {
    9751333        int iRow = which[i];
     
    10071365    }
    10081366  } else {
    1009     if (mode_) {
     1367    if (switchType==1) {
    10101368      for (i=0;i<number;i++) {
    10111369        int iRow = which[i];
     
    10581416                                       CoinIndexedVector * rowArray2)
    10591417{
     1418  if (mode_==4&&!numberSwitched_)
     1419    return;
    10601420  model_->unpack(rowArray1,sequence);
    10611421  model_->factorization()->updateColumn(rowArray2,rowArray1);
     
    10681428  int i;
    10691429
    1070   if (mode_) {
     1430  if (mode_==1) {
    10711431    for (i=0;i<number;i++) {
    10721432      int iRow = which[i];
     
    11061466  int number = numberRows + numberColumns;
    11071467  int iSequence;
    1108   if (!mode_) {
     1468  if (mode_!=1) {
    11091469    // initialize to 1.0
    11101470    // and set reference framework
  • branches/pre/ClpSimplex.cpp

    r208 r210  
    397397    }
    398398  }
    399 
    400399  arrayVector.setNumElements(number);
    401400
     
    11521151    setStatus(sequenceIn_,basic);
    11531152    if (upper_[sequenceOut_]-lower_[sequenceOut_]>0) {
    1154       if (directionOut_>0) {
     1153      // As Nonlinear costs may have moved bounds (to more feasible)
     1154      // Redo using value
     1155      if (fabs(valueOut_-lower_[sequenceOut_])<fabs(valueOut_-upper_[sequenceOut_])) {
    11551156        // going to lower
    11561157        setStatus(sequenceOut_,atLowerBound);
     
    11661167  } else {
    11671168    // flip from bound to bound
    1168     if (directionIn_==-1) {
     1169    // As Nonlinear costs may have moved bounds (to more feasible)
     1170    // Redo using value
     1171    if (fabs(valueIn_-lower_[sequenceIn_])<fabs(valueIn_-upper_[sequenceIn_])) {
    11691172      // as if from upper bound
    11701173      setStatus(sequenceIn_, atLowerBound);
     
    26632666int ClpSimplex::primal (int ifValuesPass )
    26642667{
    2665   return ((ClpSimplexPrimal *) this)->primal(ifValuesPass);
     2668  int returnCode = ((ClpSimplexPrimal *) this)->primal(ifValuesPass);
     2669  if (problemStatus_==10) {
     2670    //printf("Cleaning up with dual\n");
     2671    int savePerturbation = perturbation_;
     2672    perturbation_=100;
     2673    returnCode = ((ClpSimplexDual *) this)->dual(0);
     2674    perturbation_=savePerturbation;
     2675  }
     2676  return returnCode;
    26662677}
    26672678#include "ClpSimplexPrimalQuadratic.hpp"
     
    41694180  sequenceIn_=-1;
    41704181  sequenceOut_=-1;
     4182  secondaryStatus_=0;
    41714183
    41724184  primalTolerance_=dblParam_[ClpPrimalTolerance];
     
    41904202    factorization_->slackValue(-1.0);
    41914203    factorization_->zeroTolerance(1.0e-13);
    4192     // Switch off dense
     4204    // Switch off dense (unless special option set)
    41934205    int saveThreshold = factorization_->denseThreshold();
    41944206    factorization_->setDenseThreshold(0);
    4195 #if 0
    4196     // do perturbation if asked for
    4197 
    4198     if (perturbation_<100) {
    4199       if (algorithm_>0) {
    4200         ((ClpSimplexPrimal *) this)->perturb();
    4201       } else if (algorithm_<0) {
     4207    // If values pass then perturb (otherwise may be optimal so leave a bit)
     4208    if (ifValuesPass) {
     4209      // do perturbation if asked for
     4210     
     4211      if (perturbation_<100) {
     4212        if (algorithm_>0) {
     4213          ((ClpSimplexPrimal *) this)->perturb(0);
     4214        } else if (algorithm_<0) {
    42024215        ((ClpSimplexDual *) this)->perturb();
    4203       }
    4204     }
    4205 #endif
     4216        }
     4217      }
     4218    }
    42064219    // for primal we will change bounds using infeasibilityCost_
    42074220    if (nonLinearCost_==NULL&&algorithm_>0) {
     
    46854698  return matched;
    46864699}
     4700// Allow initial dense factorization
     4701void
     4702ClpSimplex::setInitialDenseFactorization(bool onOff)
     4703{
     4704  if (onOff)
     4705    specialOptions_ |= 8;
     4706  else
     4707    specialOptions_ &= ~8;
     4708}
     4709bool
     4710ClpSimplex::initialDenseFactorization() const
     4711{
     4712  return (specialOptions_&8)!=0;
     4713}
  • branches/pre/ClpSimplexDual.cpp

    r208 r210  
    310310      matrix_->refresh(this);
    311311      // If getting nowhere - why not give it a kick
    312 #if 0
    313312      // does not seem to work too well - do some more work
    314       if (perturbation_<101&&numberIterations_>2*(numberRows_+numberColumns_))
     313      if (perturbation_<101&&numberIterations_>2*(numberRows_+numberColumns_)) {
    315314        perturb();
    316 #endif
     315        // Can't get here if values pass
     316        gutsOfSolution(NULL,NULL);
     317        if (handler_->logLevel()>2) {
     318          handler_->message(CLP_SIMPLEX_STATUS,messages_)
     319            <<numberIterations_<<objectiveValue();
     320          handler_->printing(sumPrimalInfeasibilities_>0.0)
     321            <<sumPrimalInfeasibilities_<<numberPrimalInfeasibilities_;
     322          handler_->printing(sumDualInfeasibilities_>0.0)
     323            <<sumDualInfeasibilities_<<numberDualInfeasibilities_;
     324          handler_->printing(numberDualInfeasibilitiesWithoutFree_
     325                             <numberDualInfeasibilities_)
     326                               <<numberDualInfeasibilitiesWithoutFree_;
     327          handler_->message()<<CoinMessageEol;
     328        }
     329      }
    317330      // may factorize, checks if problem finished
    318331      statusOfProblemInDual(lastCleaned,factorType,progress_,saveDuals);
     
    23922405  }
    23932406  progressFlag_ = 0; //reset progress flag
    2394 #ifdef CLP_DEBUG
    2395   if (!rowScale_&&(handler_->logLevel()&32)) {
    2396     double * objectiveSimplex
    2397       = ClpCopyOfArray(objective(),numberColumns_,0.0);
    2398     double * rowObjectiveSimplex
    2399       = ClpCopyOfArray(rowObjective_,numberRows_,0.0);
    2400     int i;
    2401     double largest;
    2402     largest=0.0;
    2403     // direction is actually scale out not scale in
    2404     if (direction)
    2405       direction = 1.0/direction;
    2406     for (i=0;i<numberRows_;i++) {
    2407       rowObjectiveSimplex[i] *= direction;
    2408       double difference = fabs(rowObjectiveWork_[i]-rowObjectiveSimplex[i]);
    2409       if (difference>largest)
    2410         largest=difference;
    2411     }
    2412     for (i=0;i<numberColumns_;i++) {
    2413       objectiveSimplex[i] *= direction;
    2414       double difference = fabs(objectiveWork_[i]-objectiveSimplex[i]);
    2415       if (difference>largest)
    2416         largest=difference;
    2417     }
    2418     if ((handler_->logLevel()&16))
    2419       printf("difference in obj %g\n",largest);
    2420     delete [] objectiveSimplex;
    2421     delete [] rowObjectiveSimplex;
    2422   }
    2423 #endif
    24242407  handler_->message(CLP_SIMPLEX_STATUS,messages_)
    24252408    <<numberIterations_<<objectiveValue();
     
    27752758           !numberAtFakeBound()&&!numberDualInfeasibilities_) {
    27762759    problemStatus_=1;
     2760    secondaryStatus_ = 1; // and say was on cutoff
    27772761  }
    27782762  if (problemStatus_<0&&!changeMade_) {
     
    31043088  const double m1 = 0.5;
    31053089  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
    3106     if (columnLowerWork_[iColumn]<columnUpperWork_[iColumn]) {
     3090    if (columnLowerWork_[iColumn]<columnUpperWork_[iColumn]&&getStatus(iColumn)!=basic) {
    31073091      double value = perturbation;
    31083092      double currentValue = objectiveWork_[iColumn];
     
    35943578           numberAtFakeBound()) {
    35953579          returnCode=1;
     3580          secondaryStatus_ = 1; // and say was on cutoff
    35963581          // can't say anything interesting - might as well return
    35973582#ifdef CLP_DEBUG
  • branches/pre/ClpSimplexPrimal.cpp

    r207 r210  
    197197    // This says whether to restore things etc
    198198    int factorType=0;
     199    if (problemStatus_<0&&perturbation_<100) {
     200      perturb(0);
     201      // Can't get here if values pass
     202      assert (!ifValuesPass);
     203      gutsOfSolution(NULL,NULL);
     204      if (handler_->logLevel()>2) {
     205        handler_->message(CLP_SIMPLEX_STATUS,messages_)
     206          <<numberIterations_<<objectiveValue();
     207        handler_->printing(sumPrimalInfeasibilities_>0.0)
     208          <<sumPrimalInfeasibilities_<<numberPrimalInfeasibilities_;
     209        handler_->printing(sumDualInfeasibilities_>0.0)
     210          <<sumDualInfeasibilities_<<numberDualInfeasibilities_;
     211        handler_->printing(numberDualInfeasibilitiesWithoutFree_
     212                           <numberDualInfeasibilities_)
     213                             <<numberDualInfeasibilitiesWithoutFree_;
     214        handler_->message()<<CoinMessageEol;
     215      }
     216    }
    199217    /*
    200218      Status of problem:
     
    223241      matrix_->refresh(this);
    224242      // If getting nowhere - why not give it a kick
    225 #if 0
    226       // primal perturbation not coded yet
     243#if 1
    227244      if (perturbation_<101&&numberIterations_>2*(numberRows_+numberColumns_))
    228         perturb();
     245        perturb(1);
    229246#endif
    230247      // If we have done no iterations - special
     
    487504    gutsOfSolution(NULL,NULL);
    488505  }
     506  // Flag to say whether to go to dual to clean up
     507  bool goToDual=false;
    489508  // really for free variables in
    490509  //if((progressFlag_&2)!=0)
     
    502521                       <<numberDualInfeasibilitiesWithoutFree_;
    503522  handler_->message()<<CoinMessageEol;
    504   assert (primalFeasible());
     523  if (!primalFeasible()) {
     524    nonLinearCost_->checkInfeasibilities(primalTolerance_);
     525    gutsOfSolution(NULL,NULL);
     526  }
    505527  // we may wish to say it is optimal even if infeasible
    506528  bool alwaysOptimal = (specialOptions_&1)!=0;
     
    550572        if ((infeasibilityCost_>=1.0e18||
    551573             numberDualInfeasibilities_==0)&&perturbation_==101) {
    552           unPerturb(); // stop any further perturbation
     574          goToDual=unPerturb(); // stop any further perturbation
    553575          numberDualInfeasibilities_=1; // carry on
    554576        }
     
    590612      // may be optimal
    591613      if (perturbation_==101) {
    592         unPerturb(); // stop any further perturbation
     614        goToDual=unPerturb(); // stop any further perturbation
    593615        lastCleaned=-1; // carry on
    594616      }
    595       if ( lastCleaned!=numberIterations_||unflag()) {
     617      bool unflagged = unflag();
     618      if ( lastCleaned!=numberIterations_||unflagged) {
    596619        handler_->message(CLP_PRIMAL_OPTIMAL,messages_)
    597620          <<primalTolerance_
     
    651674#endif
    652675          gutsOfSolution(NULL,NULL);
    653           problemStatus_ = -1;
     676          if (sumOfRelaxedDualInfeasibilities_ == 0.0&&
     677              sumOfRelaxedPrimalInfeasibilities_ == 0.0) {
     678            // say optimal (with these bounds etc)
     679            numberDualInfeasibilities_ = 0;
     680            sumDualInfeasibilities_ = 0.0;
     681            numberPrimalInfeasibilities_ = 0;
     682            sumPrimalInfeasibilities_ = 0.0;
     683          }
     684          if (dualFeasible()&&!nonLinearCost_->numberInfeasibilities()&&lastCleaned>=0)
     685            problemStatus_=0;
     686          else
     687            problemStatus_ = -1;
    654688        } else {
    655689          problemStatus_=0; // optimal
     
    670704          // back off weight
    671705          infeasibilityCost_ = 1.0e13;
    672           unPerturb(); // stop any further perturbation
     706          goToDual=unPerturb(); // stop any further perturbation
    673707        }
    674708        //we need infeasiblity cost changed
     
    726760  }
    727761  lastGoodIteration_ = numberIterations_;
     762  if (goToDual)
     763    problemStatus_=10; // try dual
    728764#if 0
    729765  double thisObj = progress->lastObjective(0);
     
    861897  int iIndex;
    862898#if 0
    863   if (numberIterations_==1147)
     899  if (numberIterations_<=39)
     900    handler_->setLogLevel(63);
     901  else
     902    handler_->setLogLevel(2);
     903  if (numberIterations_==38)
    864904    printf("trouble\n");
     905  assert (solution_[29176]>-1.0e20);
    865906#endif
    866907  for (iIndex=0;iIndex<number;iIndex++) {
     
    13241365// Perturbs problem
    13251366void
    1326 ClpSimplexPrimal::perturb()
     1367ClpSimplexPrimal::perturb(int type)
    13271368{
    13281369  if (perturbation_>100)
     
    13321373  double perturbation=1.0e-20;
    13331374  // maximum fraction of rhs/bounds to perturb
    1334   double maximumFraction = 1.0e-4;
     1375  double maximumFraction = 1.0e-8;
    13351376  if (perturbation_>=50) {
    13361377    perturbation = 1.0e-4;
     
    13511392      }
    13521393    }
    1353     perturbation *= 1.0e-8;
    13541394  } else if (perturbation_<100) {
    13551395    perturbation = pow(10.0,perturbation_);
    13561396    // user is in charge
    1357     maximumFraction = 1.0e100;
    1358   }
     1397    maximumFraction = 1.0;
     1398  }
     1399  double largestZero=0.0;
     1400  double largest=0.0;
     1401  double largestPerCent=0.0;
     1402  bool printOut=(handler_->logLevel()==63);
     1403  // Check if all slack
     1404  int number=0;
     1405  int iSequence;
     1406  for (iSequence=0;iSequence<numberRows_;iSequence++) {
     1407    if (getRowStatus(iSequence)==basic)
     1408      number++;
     1409  }
     1410  if (number!=numberRows_)
     1411    type=1;
    13591412  // modify bounds
     1413  if (type==1) {
     1414    double multiplier = perturbation*maximumFraction;
     1415    for (iSequence=0;iSequence<numberRows_+numberColumns_;iSequence++) {
     1416      if (getStatus(iSequence)==basic) {
     1417        double solutionValue = solution_[iSequence];
     1418        double lowerValue = lower_[iSequence];
     1419        double upperValue = upper_[iSequence];
     1420        double difference = upperValue-lowerValue;
     1421        difference = min(difference,perturbation);
     1422        difference = min(difference,fabs(solutionValue)+1.0);
     1423        double value = CoinDrand48()*multiplier*(difference+1.0);
     1424        if (solutionValue-lowerValue<=primalTolerance_) {
     1425          lower_[iSequence] -= value;
     1426        } else if (upperValue-solutionValue<=primalTolerance_) {
     1427          upper_[iSequence] += value;
     1428        } else {
     1429#if 0
     1430          if (iSequence>=numberColumns_) {
     1431            // may not be at bound - but still perturb (unless free)
     1432            if (upperValue>1.0e30&&lowerValue<-1.0e30)
     1433              value=0.0;
     1434            else
     1435              value = - value; // as -1.0 in matrix
     1436          } else {
     1437            value = 0.0;
     1438          }
     1439#else
     1440          value=0.0;
     1441#endif
     1442        }
     1443        if (value) {
     1444          if (printOut)
     1445            printf("col %d lower from %g to %g, upper from %g to %g\n",
     1446                   iSequence,lower_[iSequence],lowerValue,upper_[iSequence],upperValue);
     1447          if (solutionValue) {
     1448            largest = max(largest,value);
     1449            if (value>(fabs(solutionValue)+1.0)*largestPerCent)
     1450              largestPerCent=value/(fabs(solutionValue)+1.0);
     1451          } else {
     1452            largestZero = max(largestZero,value);
     1453          }
     1454        }
     1455      }
     1456    }
     1457  } else {
     1458    for (i=0;i<numberColumns_;i++) {
     1459      double lowerValue=lower_[i], upperValue=upper_[i];
     1460      if (upperValue>lowerValue+primalTolerance_) {
     1461        double value = CoinDrand48()*perturbation*maximumFraction;
     1462        if (lowerValue>-1.0e20&&lowerValue)
     1463          lowerValue -= value * (max(1.0,1.0e-5*fabs(lowerValue)));
     1464        if (upperValue<1.0e20&&upperValue)
     1465          upperValue += value * (max(1.0,1.0e-5*fabs(upperValue)));
     1466        if (lowerValue!=lower_[i]) {
     1467          double difference = fabs(lowerValue-lower_[i]);
     1468          largest = max(largest,difference);
     1469          if (difference>fabs(lower_[i])*largestPerCent)
     1470            largestPerCent=fabs(difference/lower_[i]);
     1471        }
     1472        if (upperValue!=upper_[i]) {
     1473          double difference = fabs(upperValue-upper_[i]);
     1474          largest = max(largest,difference);
     1475          if (difference>fabs(upper_[i])*largestPerCent)
     1476            largestPerCent=fabs(difference/upper_[i]);
     1477        }
     1478        if (printOut)
     1479          printf("col %d lower from %g to %g, upper from %g to %g\n",
     1480                 i,lower_[i],lowerValue,upper_[i],upperValue);
     1481      }
     1482      if (solution_[i]==lower_[i])
     1483        solution_[i]=lowerValue;
     1484      else if (solution_[i]==upper_[i])
     1485        solution_[i]=upperValue;
     1486      lower_[i]=lowerValue;
     1487      upper_[i]=upperValue;
     1488    }
     1489    for (;i<numberColumns_+numberRows_;i++) {
     1490      double lowerValue=lower_[i], upperValue=upper_[i];
     1491      double value = CoinDrand48()*perturbation*maximumFraction;
     1492      if (upperValue>lowerValue+primalTolerance_) {
     1493        if (lowerValue>-1.0e20&&lowerValue)
     1494          lowerValue -= value * (max(1.0,1.0e-5*fabs(lowerValue)));
     1495        if (upperValue<1.0e20&&upperValue)
     1496          upperValue += value * (max(1.0,1.0e-5*fabs(upperValue)));
     1497      } else if (upperValue>0.0) {
     1498        upperValue -= value * (max(1.0,1.0e-5*fabs(lowerValue)));
     1499        lowerValue -= value * (max(1.0,1.0e-5*fabs(lowerValue)));
     1500      } else if (upperValue<0.0) {
     1501        upperValue += value * (max(1.0,1.0e-5*fabs(lowerValue)));
     1502        lowerValue += value * (max(1.0,1.0e-5*fabs(lowerValue)));
     1503      } else {
     1504      }
     1505      if (lowerValue!=lower_[i]) {
     1506        double difference = fabs(lowerValue-lower_[i]);
     1507        largest = max(largest,difference);
     1508        if (difference>fabs(lower_[i])*largestPerCent)
     1509          largestPerCent=fabs(difference/lower_[i]);
     1510      }
     1511      if (upperValue!=upper_[i]) {
     1512        double difference = fabs(upperValue-upper_[i]);
     1513        largest = max(largest,difference);
     1514        if (difference>fabs(upper_[i])*largestPerCent)
     1515          largestPerCent=fabs(difference/upper_[i]);
     1516      }
     1517      if (printOut)
     1518        printf("row %d lower from %g to %g, upper from %g to %g\n",
     1519               i-numberColumns_,lower_[i],lowerValue,upper_[i],upperValue);
     1520      if (solution_[i]==lower_[i])
     1521        solution_[i]=lowerValue;
     1522      else if (solution_[i]==upper_[i])
     1523        solution_[i]=upperValue;
     1524      lower_[i]=lowerValue;
     1525      upper_[i]=upperValue;
     1526    }
     1527  }
    13601528  handler_->message(CLP_SIMPLEX_PERTURB,messages_)
    1361     <<perturbation
     1529    <<100.0*maximumFraction<<perturbation<<largest<<100.0*largestPerCent<<largestZero
    13621530    <<CoinMessageEol;
    1363   for (i=0;i<numberColumns_;i++) {
    1364     double lowerValue=lower_[i], upperValue=upper_[i];
    1365     if (upperValue>lowerValue+primalTolerance_) {
    1366       double value = CoinDrand48()*perturbation;
    1367       if (lowerValue>-1.0e20&&lowerValue)
    1368         lowerValue -= value * (max(1.0,1.0e-5*fabs(lowerValue)));
    1369       if (upperValue<1.0e20&&upperValue)
    1370         upperValue += value * (max(1.0,1.0e-5*fabs(upperValue)));
    1371     }
    1372     lower_[i]=lowerValue;
    1373     upper_[i]=upperValue;
    1374   }
    1375   for (;i<numberColumns_+numberRows_;i++) {
    1376     double lowerValue=lower_[i], upperValue=upper_[i];
    1377     double value = CoinDrand48()*perturbation;
    1378     if (upperValue>lowerValue+primalTolerance_) {
    1379       if (lowerValue>-1.0e20)
    1380         lowerValue -= value * (max(1.0,1.0e-5*fabs(lowerValue)));
    1381       if (upperValue<1.0e20)
    1382         upperValue += value * (max(1.0,1.0e-5*fabs(upperValue)));
    1383     } else if (upperValue>0.0) {
    1384         lowerValue -= value * (max(1.0,1.0e-5*fabs(lowerValue)));
    1385         upperValue -= value * (max(1.0,1.0e-5*fabs(lowerValue)));
    1386     } else if (upperValue<0.0) {
    1387         lowerValue += value * (max(1.0,1.0e-5*fabs(lowerValue)));
    1388         upperValue += value * (max(1.0,1.0e-5*fabs(lowerValue)));
    1389     } else {
    1390     }
    1391     lower_[i]=lowerValue;
    1392     upper_[i]=upperValue;
    1393   }
     1531  // redo nonlinear costs
    13941532  // say perturbed
    13951533  perturbation_=101;
    13961534}
    13971535// un perturb
    1398 void
     1536bool
    13991537ClpSimplexPrimal::unPerturb()
    14001538{
     1539  if (perturbation_!=101)
     1540    return false;
    14011541  // put back original bounds and costs
    14021542  createRim(7);
     
    14091549  // move non basic variables to new bounds
    14101550  nonLinearCost_->checkInfeasibilities(0.0);
     1551#if 1
     1552  // Try using dual
     1553  return true;
     1554#else
    14111555  gutsOfSolution(NULL,NULL);
     1556  return false;
     1557#endif
     1558 
    14121559}
    14131560// Unflag all variables and return number unflagged
     
    17291876        valueOut_ = upperOut_;
    17301877      }
    1731       double lowerValue = lower_[sequenceOut_];
    1732       double upperValue = upper_[sequenceOut_];
    1733       assert(valueOut_>=lowerValue-primalTolerance_&&
    1734              valueOut_<=upperValue+primalTolerance_);
     1878      assert(valueOut_>=lower_[sequenceOut_]-primalTolerance_&&
     1879             valueOut_<=upper_[sequenceOut_]+primalTolerance_);
    17351880      // may not be exactly at bound and bounds may have changed
    17361881      // Make sure outgoing looks feasible
    1737       double useValue=valueOut_;
    1738       if (valueOut_<=lowerValue+primalTolerance_) {
    1739         directionOut_=1;
    1740         useValue= lower_[sequenceOut_];
    1741       } else if (valueOut_>=upperValue-primalTolerance_) {
    1742         directionOut_=-1;
    1743         useValue= upper_[sequenceOut_];
    1744       } else {
    1745         printf("*** variable wandered off bound %g %g %g!\n",
    1746                lowerValue,valueOut_,upperValue);
    1747         if (valueOut_-lowerValue<=upperValue-valueOut_) {
    1748           useValue= lower_[sequenceOut_];
    1749           valueOut_=lowerValue;
    1750           directionOut_=1;
    1751         } else {
    1752           useValue= upper_[sequenceOut_];
    1753           valueOut_=upperValue;
    1754           directionOut_=-1;
    1755         }
    1756       }
     1882      directionOut_=nonLinearCost_->setOneOutgoing(sequenceOut_,valueOut_);
    17571883      solution_[sequenceOut_]=valueOut_;
    1758       nonLinearCost_->setOne(sequenceOut_,useValue);
    17591884    }
    17601885    // change cost and bounds on incoming if primal
  • branches/pre/ClpSimplexPrimalQuadratic.cpp

    r205 r210  
    11911191                 valueOut_<=upperValue+primalTolerance_);
    11921192          // may not be exactly at bound and bounds may have changed
     1193#if 1
    11931194          // Make sure outgoing looks feasible
    11941195          double useValue=valueOut_;
     
    12141215          solution_[sequenceOut_]=valueOut_;
    12151216          nonLinearCost_->setOne(sequenceOut_,useValue);
     1217#else
     1218          directionOut_=nonLinearCost_->setOneOutgoing(sequenceOut_,valueOut_);
     1219          solution_[sequenceOut_]=valueOut_;
     1220#endif
    12161221        }
    12171222        // change cost and bounds on incoming if primal
  • branches/pre/ClpSolve.cpp

    r208 r210  
    1010
    1111#include "CoinHelperFunctions.hpp"
     12#include "CoinSort.hpp"
    1213#include "ClpFactorization.hpp"
    1314#include "ClpSimplex.hpp"
     15#include "ClpSolve.hpp"
    1416#include "ClpPackedMatrix.hpp"
    1517#include "ClpPlusMinusOneMatrix.hpp"
    1618#include "ClpMessage.hpp"
     19#include "CoinTime.hpp"
    1720
    1821#include "ClpPresolve.hpp"
     
    4245 */
    4346int
    44 ClpSimplex::initialSolve(SolveType method, PresolveType presolve,
    45                          int specialOption)
    46 {
     47ClpSimplex::initialSolve(ClpSolve & options)
     48{
     49  ClpSolve::SolveType method=options.getSolveType();
     50  ClpSolve::PresolveType presolve = options.getPresolveType();
    4751  int saveMaxIterations = maximumIterations();
    4852  int finalStatus=-1;
     
    5357  ClpMatrixBase * saveMatrix=NULL;
    5458  ClpSimplex * model2 = this;
    55   bool interrupt = ((specialOption&16)==0);
     59  bool interrupt = (options.getSpecialOption(2)==0);
    5660  sighandler_t saveSignal=SIG_DFL;
    5761  if (interrupt) {
     
    6569  double timeCore=0.0;
    6670  double timeSimplex=0.0;
    67   assert (method!=automatic); // later
    6871  int savePerturbation=perturbation_;
    69   if ((specialOption&1)!=0||method==usePrimal)
    70     perturbation_=100;
    7172  int saveScaling = scalingFlag_;
    72   if ((specialOption&2)!=0)
    73     scalingFlag_=0;
    74 
    75   if (presolve!=presolveOff) {
     73
     74  if (presolve!=ClpSolve::presolveOff) {
    7675    int numberPasses=5;
    77     if (presolve==presolveMaximum)
    78       numberPasses=25;
     76    if (presolve==ClpSolve::presolveNumber) {
     77      numberPasses=options.getPresolvePasses();
     78      presolve = ClpSolve::presolveOn;
     79    }
    7980    model2 = pinfo.presolvedModel(*this,1.0e-8,
    8081                                  false,numberPasses,true);
     
    8687    timeX=time2;
    8788    if (model2) {
    88       if (method==useDual) {
     89      if (method==ClpSolve::useDual) {
    8990        int numberInfeasibilities = model2->tightenPrimalBounds();
    9091        if (numberInfeasibilities) {
     
    9293            <<CoinMessageEol;
    9394          model2 = this;
    94           presolve=presolveOff;
     95          presolve=ClpSolve::presolveOff;
    9596        }
    9697      }
     
    99100        <<CoinMessageEol;
    100101      model2 = this;
    101       presolve=presolveOff;
     102      presolve=ClpSolve::presolveOff;
     103    }
     104    // We may be better off using original
     105    if (numberRows_<1.01*model2->numberRows_&&numberColumns_<1.01*model2->numberColumns_) {
     106      delete model2;
     107      model2 = this;
     108      presolve=ClpSolve::presolveOff;
    102109    }
    103110  }
     
    106113  // See if worth trying +- one matrix
    107114  bool plusMinus=false;
    108   if ((specialOption&32)==0) {
    109     int numberElements=model2->getNumElements();
     115  int numberElements=model2->getNumElements();
     116  // For below >0 overrides
     117  // 0 means no, -1 means maybe
     118  int doIdiot=0;
     119  int doCrash=0;
     120  int doSprint=0;
     121  switch (options.getSpecialOption(1)) {
     122  case 0:
     123    doIdiot=-1;
     124    doCrash=-1;
     125    doSprint=-1;
     126    break;
     127  case 1:
     128    doIdiot=0;
     129    doCrash=1;
     130    doSprint=0;
     131    break;
     132  case 2:
     133    doIdiot=1;
     134    doCrash=0;
     135    doSprint=0;
     136    break;
     137  case 3:
     138    doIdiot=0;
     139    doCrash=0;
     140    doSprint=1;
     141    break;
     142  case 4:
     143    doIdiot=0;
     144    doCrash=0;
     145    doSprint=0;
     146    break;
     147  case 5:
     148    doIdiot=0;
     149    doCrash=-1;
     150    doSprint=-1;
     151    break;
     152  case 6:
     153    doIdiot=-1;
     154    doCrash=-1;
     155    doSprint=0;
     156    break;
     157  case 7:
     158    doIdiot=-1;
     159    doCrash=0;
     160    doSprint=-1;
     161    break;
     162  case 8:
     163    doIdiot=-1;
     164    doCrash=0;
     165    doSprint=0;
     166    break;
     167  case 9:
     168    doIdiot=0;
     169    doCrash=0;
     170    doSprint=-1;
     171    break;
     172  default:
     173    abort();
     174  }
     175  int numberColumns = model2->numberColumns();
     176  int numberRows = model2->numberRows();
     177  // If not all slack basis - switch off all
     178  int number=0;
     179  int iRow;
     180  for (iRow=0;iRow<numberRows;iRow++)
     181    if (model2->getRowStatus(iRow)==basic)
     182      number++;
     183  if (number<numberRows) {
     184    doIdiot=0;
     185    doCrash=0;
     186    doSprint=0;
     187  }
     188  if (options.getSpecialOption(3)==0) {
    110189    if(numberElements>100000)
    111190      plusMinus=true;
    112     if(numberElements>10000&&(specialOption&12)==0&&method==usePrimal)
     191    if(numberElements>10000&&((doIdiot||doSprint)&&method==ClpSolve::usePrimal))
    113192      plusMinus=true;
    114193  }
     
    138217    model2->setFactorizationFrequency(100+model2->numberRows()/100);
    139218  }
    140   int numberColumns = model2->numberColumns();
    141   int numberRows = model2->numberRows();
    142   if (method==useDual) {
    143     if ((specialOption&4)!=0)
     219  if (method==ClpSolve::usePrimalorSprint) {
     220    if (doSprint<0) {
     221      if(numberRows*10>numberColumns||numberColumns<6000
     222         ||(numberRows*20>numberColumns&&!plusMinus))
     223        method=ClpSolve::usePrimal; // switch off sprint
     224    } else if (doSprint==0) {
     225      method=ClpSolve::usePrimal; // switch off sprint
     226    }
     227  }
     228  if (method==ClpSolve::useDual) {
     229    if (options.getSpecialOption(0)!=0)
    144230      model2->crash(1000,1);
    145231    model2->dual();
     
    150236      <<CoinMessageEol;
    151237    timeX=time2;
    152   } else {
     238  } else if (method==ClpSolve::usePrimal) {
    153239#ifdef CLP_IDIOT
    154     if ((specialOption&12)==0) {
     240    if (doIdiot) {
    155241      int nPasses=0;
    156       if (numberRows>2000&&numberColumns>2*numberRows) {
     242      Idiot info(*model2);
     243      if (numberRows>1000&&numberColumns>2*numberRows) {
    157244        if (plusMinus) {
    158           nPasses = 10+numberColumns/1000;
    159           nPasses = min(nPasses,100);
     245          // look at rhs
     246          int iRow;
     247          double largest=0.0;
     248          double smallest = 1.0e30;
     249          for (iRow=0;iRow<numberRows;iRow++) {
     250            double value;
     251            value = model2->rowLower_[iRow];
     252            if (value) {
     253              largest = max(largest,fabs(value));
     254              smallest=min(smallest,fabs(value));
     255            }
     256            value = model2->rowUpper_[iRow];
     257            if (value) {
     258              largest = max(largest,fabs(value));
     259              smallest=min(smallest,fabs(value));
     260            }
     261          }
     262          if (largest/smallest>2.0) {
     263            nPasses = 10+numberColumns/100000;
     264            nPasses = min(nPasses,50);
     265            nPasses = max(nPasses,15);
     266            if (numberElements<3*numberColumns)
     267              nPasses=0; // probably not worh it
     268          } else if (largest/smallest>1.01||numberElements<=3*numberColumns) {
     269            nPasses = 10+numberColumns/1000;
     270            nPasses = min(nPasses,100);
     271            nPasses = max(nPasses,30);
     272          } else {
     273            nPasses = 10+numberColumns/1000;
     274            nPasses = min(nPasses,200);
     275            nPasses = max(nPasses,100);
     276            info.setStartingWeight(1.0e-1);
     277            info.setReduceIterations(6);
     278            //info.setFeasibilityTolerance(1.0e-7);
     279          }
    160280        } else {
    161281          nPasses = 10+numberColumns/100000;
     
    164284          else
    165285            nPasses=5;
     286          if (numberElements<3*numberColumns)
     287            nPasses=0; // probably not worh it
     288          else
     289            nPasses = max(nPasses,5);
     290          //info.setStartingWeight(1.0e-1);
    166291        }
    167292      }
     293      if (doIdiot>0) {
     294        // pick up number passes
     295        nPasses=options.getExtraInfo(1);
     296        if (nPasses>70)
     297          info.setReduceIterations(6);
     298      }
    168299      if (nPasses) {
    169         Idiot info(*model2);
     300        doCrash=0;
    170301        info.crash(nPasses,model2->messageHandler(),model2->messagesPointer());
    171302        time2 = CoinCpuTime();
     
    178309    }
    179310#endif
    180     if ((specialOption&4)!=0)
     311    // ?
     312    if (doCrash)
    181313      model2->crash(1000,1);
    182314    model2->primal(1);
     
    188320      <<CoinMessageEol;
    189321    timeX=time2;
     322  } else if (method==ClpSolve::usePrimalorSprint) {
     323    // Sprint
     324    /*
     325      This driver implements what I called Sprint when I introduced the idea
     326      many years ago.  Cplex calls it "sifting" which I think is just as silly.
     327      When I thought of this trivial idea
     328      it reminded me of an LP code of the 60's called sprint which after
     329      every factorization took a subset of the matrix into memory (all
     330      64K words!) and then iterated very fast on that subset.  On the
     331      problems of those days it did not work very well, but it worked very
     332      well on aircrew scheduling problems where there were very large numbers
     333      of columns all with the same flavor.
     334    */
     335   
     336    /* The idea works best if you can get feasible easily.  To make it
     337       more general we can add in costed slacks */
     338   
     339    int originalNumberColumns = model2->numberColumns();
     340    int numberRows = model2->numberRows();
     341   
     342    // We will need arrays to choose variables.  These are too big but ..
     343    double * weight = new double [numberRows+originalNumberColumns];
     344    int * sort = new int [numberRows+originalNumberColumns];
     345    int numberSort=0;
     346    // We are going to add slacks to get feasible.
     347    // initial list will just be artificials
     348    // first we will set all variables as close to zero as possible
     349    int iColumn;
     350    const double * columnLower = model2->columnLower();
     351    const double * columnUpper = model2->columnUpper();
     352    double * columnSolution = model2->primalColumnSolution();
     353   
     354    for (iColumn=0;iColumn<originalNumberColumns;iColumn++) {
     355      double value =0.0;
     356      if (columnLower[iColumn]>0.0)
     357        value = columnLower[iColumn];
     358      else if (columnUpper[iColumn]<0.0)
     359        value = columnUpper[iColumn];
     360      columnSolution[iColumn]=value;
     361    }
     362    // now see what that does to row solution
     363    double * rowSolution = model2->primalRowSolution();
     364    memset (rowSolution,0,numberRows*sizeof(double));
     365    model2->times(1.0,columnSolution,rowSolution);
     366   
     367    int * addStarts = new int [numberRows+1];
     368    int * addRow = new int[numberRows];
     369    double * addElement = new double[numberRows];
     370    const double * lower = model2->rowLower();
     371    const double * upper = model2->rowUpper();
     372    addStarts[0]=0;
     373    int numberArtificials=0;
     374    double * addCost = new double [numberRows];
     375    const double penalty=1.0e8;
     376    int iRow;
     377    for (iRow=0;iRow<numberRows;iRow++) {
     378      if (lower[iRow]>rowSolution[iRow]) {
     379        addRow[numberArtificials]=iRow;
     380        addElement[numberArtificials]=1.0;
     381        addCost[numberArtificials]=penalty;
     382        numberArtificials++;
     383        addStarts[numberArtificials]=numberArtificials;
     384      } else if (upper[iRow]<rowSolution[iRow]) {
     385        addRow[numberArtificials]=iRow;
     386        addElement[numberArtificials]=-1.0;
     387        addCost[numberArtificials]=penalty;
     388        numberArtificials++;
     389        addStarts[numberArtificials]=numberArtificials;
     390      }
     391    }
     392    model2->addColumns(numberArtificials,NULL,NULL,addCost,
     393                       addStarts,addRow,addElement);
     394    delete [] addStarts;
     395    delete [] addRow;
     396    delete [] addElement;
     397    delete [] addCost;
     398    // look at rhs to see if to perturb
     399    double largest=0.0;
     400    double smallest = 1.0e30;
     401    for (iRow=0;iRow<numberRows;iRow++) {
     402      double value;
     403      value = fabs(model2->rowLower_[iRow]);
     404      if (value&&value<1.0e30) {
     405        largest = max(largest,value);
     406        smallest=min(smallest,value);
     407      }
     408      value = fabs(model2->rowUpper_[iRow]);
     409      if (value&&value<1.0e30) {
     410        largest = max(largest,value);
     411        smallest=min(smallest,value);
     412      }
     413    }
     414    double * saveLower = NULL;
     415    double * saveUpper = NULL;
     416    if (largest<2.01*smallest) {
     417      // perturb - so switch off standard
     418      model2->setPerturbation(100);
     419      saveLower = new double[numberRows];
     420      memcpy(saveLower,model2->rowLower_,numberRows*sizeof(double));
     421      saveUpper = new double[numberRows];
     422      memcpy(saveUpper,model2->rowUpper_,numberRows*sizeof(double));
     423      double * lower = model2->rowLower();
     424      double * upper = model2->rowUpper();
     425      for (iRow=0;iRow<numberRows;iRow++) {
     426        double lowerValue=lower[iRow], upperValue=upper[iRow];
     427        double value = CoinDrand48();
     428        if (upperValue>lowerValue+primalTolerance_) {
     429          if (lowerValue>-1.0e20&&lowerValue)
     430            lowerValue -= value * 1.0e-4*fabs(lowerValue);
     431          if (upperValue<1.0e20&&upperValue)
     432            upperValue += value * 1.0e-4*fabs(upperValue);
     433        } else if (upperValue>0.0) {
     434          upperValue -= value * 1.0e-4*fabs(lowerValue);
     435          lowerValue -= value * 1.0e-4*fabs(lowerValue);
     436        } else if (upperValue<0.0) {
     437          upperValue += value * 1.0e-4*fabs(lowerValue);
     438          lowerValue += value * 1.0e-4*fabs(lowerValue);
     439        } else {
     440        }
     441        lower[iRow]=lowerValue;
     442        upper[iRow]=upperValue;
     443      }
     444    }
     445    int i;
     446    // Set up initial list
     447    if (numberArtificials) {
     448      numberSort=numberArtificials;
     449      for (i=0;i<numberSort;i++)
     450        sort[i] = i+originalNumberColumns;
     451    } else {
     452      numberSort = min(numberRows_,numberColumns_);
     453      for (i=0;i<numberSort;i++)
     454        sort[i] = i;
     455    }
     456   
     457    // redo as will have changed
     458    columnLower = model2->columnLower();
     459    columnUpper = model2->columnUpper();
     460    int numberColumns = model2->numberColumns();
     461    double * fullSolution = model2->primalColumnSolution();
     462   
     463    // Just do this number of passes
     464    int maxPass=100;
     465    if (doSprint>0)
     466      maxPass=options.getExtraInfo(1);
     467    int iPass;
     468    double lastObjective=1.0e31;
     469    // It will be safe to allow dense
     470    model2->setInitialDenseFactorization(true);
     471   
     472    // Just take this number of columns in small problem
     473    int smallNumberColumns = min(3*numberRows,numberColumns);
     474    smallNumberColumns = max(smallNumberColumns,3000);
     475    // We will be using all rows
     476    int * whichRows = new int [numberRows];
     477    for (int iRow=0;iRow<numberRows;iRow++)
     478      whichRows[iRow]=iRow;
     479    double originalOffset;
     480    model2->getDblParam(ClpObjOffset,originalOffset);
     481    for (iPass=0;iPass<maxPass;iPass++) {
     482      //printf("Bug until submodel new version\n");
     483      //CoinSort_2(sort,sort+numberSort,weight);
     484      // Create small problem
     485      ClpSimplex small(model2,numberRows,whichRows,numberSort,sort);
     486      small.setPerturbation(model2->perturbation());
     487      // now see what variables left out do to row solution
     488      double * rowSolution = model2->primalRowSolution();
     489      double * sumFixed = new double[numberRows];
     490      memset (sumFixed,0,numberRows*sizeof(double));
     491      int iRow,iColumn;
     492      // zero out ones in small problem
     493      for (iColumn=0;iColumn<numberSort;iColumn++) {
     494        int kColumn = sort[iColumn];
     495        fullSolution[kColumn]=0.0;
     496      }
     497      // Get objective offset
     498      double offset=0.0;
     499      const double * objective = model2->objective();
     500      for (iColumn=0;iColumn<numberColumns_;iColumn++)
     501        offset += fullSolution[iColumn]*objective[iColumn];
     502      small.setDblParam(ClpObjOffset,originalOffset-offset);
     503      model2->times(1.0,fullSolution,sumFixed);
     504     
     505      double * lower = small.rowLower();
     506      double * upper = small.rowUpper();
     507      for (iRow=0;iRow<numberRows;iRow++) {
     508        if (lower[iRow]>-1.0e50)
     509          lower[iRow] -= sumFixed[iRow];
     510        if (upper[iRow]<1.0e50)
     511          upper[iRow] -= sumFixed[iRow];
     512        rowSolution[iRow] -= sumFixed[iRow];
     513      }
     514      delete [] sumFixed;
     515      // Solve
     516      if (interrupt)
     517        currentModel = &small;
     518      small.primal();
     519      // move solution back
     520      const double * solution = small.primalColumnSolution();
     521      for (iColumn=0;iColumn<numberSort;iColumn++) {
     522        int kColumn = sort[iColumn];
     523        model2->setColumnStatus(kColumn,small.getColumnStatus(iColumn));
     524        fullSolution[kColumn]=solution[iColumn];
     525      }
     526      for (iRow=0;iRow<numberRows;iRow++)
     527        model2->setRowStatus(iRow,small.getRowStatus(iRow));
     528      memcpy(model2->primalRowSolution(),small.primalRowSolution(),
     529             numberRows*sizeof(double));
     530      // get reduced cost for large problem
     531      memcpy(weight,model2->objective(),numberColumns*sizeof(double));
     532      model2->transposeTimes(-1.0,small.dualRowSolution(),weight);
     533      int numberNegative=0;
     534      double sumNegative = 0.0;
     535      // now massage weight so all basic in plus good djs
     536      for (iColumn=0;iColumn<numberColumns;iColumn++) {
     537        double dj = weight[iColumn]*optimizationDirection_;
     538        double value = fullSolution[iColumn];
     539        if (model2->getColumnStatus(iColumn)==ClpSimplex::basic)
     540          dj = -1.0e50;
     541        else if (dj<0.0&&value<columnUpper[iColumn])
     542          dj = dj;
     543        else if (dj>0.0&&value>columnLower[iColumn])
     544          dj = -dj;
     545        else if (columnUpper[iColumn]>columnLower[iColumn])
     546          dj = fabs(dj);
     547        else
     548          dj = 1.0e50;
     549        weight[iColumn] = dj;
     550        if (dj<-dualTolerance_&&dj>-1.0e50) {
     551          numberNegative++;
     552          sumNegative -= dj;
     553        }
     554        sort[iColumn] = iColumn;
     555      }
     556      handler_->message(CLP_SPRINT,messages_)
     557        <<iPass+1<<small.numberIterations()<<small.objectiveValue()<<sumNegative
     558        <<numberNegative
     559        <<CoinMessageEol;
     560      if ((small.objectiveValue()>lastObjective-1.0e-7&&iPass>5)||
     561          !small.numberIterations()||
     562          iPass==maxPass-1||small.status()==3) {
     563       
     564        break; // finished
     565      } else {
     566        lastObjective = small.objectiveValue();
     567        // sort
     568        CoinSort_2(weight,weight+numberColumns,sort);
     569        numberSort = smallNumberColumns;
     570      }
     571    }
     572    if (interrupt)
     573      currentModel = model2;
     574    for (i=0;i<numberArtificials;i++)
     575      sort[i] = i + originalNumberColumns;
     576    model2->deleteColumns(numberArtificials,sort);
     577    delete [] weight;
     578    delete [] sort;
     579    delete [] whichRows;
     580    if (saveLower) {
     581      // unperturb
     582      memcpy(model2->rowLower_,saveLower,numberRows*sizeof(double));
     583      delete [] saveLower;
     584      memcpy(model2->rowUpper_,saveUpper,numberRows*sizeof(double));
     585      delete [] saveUpper;
     586    }
     587    model2->primal(1);
     588    model2->setPerturbation(savePerturbation);
     589    time2 = CoinCpuTime();
     590    timeCore = time2-timeX;
     591    timeSimplex = timeCore;
     592    handler_->message(CLP_INTERVAL_TIMING,messages_)
     593      <<"Sprint"<<timeCore<<time2-time1
     594      <<CoinMessageEol;
     595    timeX=time2;
     596  } else {
     597    assert (method!=ClpSolve::automatic); // later
     598    time2=0.0;
    190599  }
    191600  if (saveMatrix) {
     
    196605  numberIterations = model2->numberIterations();
    197606  finalStatus=model2->status();
    198   if (presolve==presolveOn) {
     607  if (presolve==ClpSolve::presolveOn) {
    199608    int saveLevel = logLevel();
    200609    setLogLevel(1);
     
    231640  handler_->message(CLP_TIMING,messages_)
    232641    <<statusMessage[finalStatus+1]<<objectiveValue()<<numberIterations<<time2-time1;
    233   handler_->printing(presolve==presolveOn)
     642  handler_->printing(presolve==ClpSolve::presolveOn)
    234643    <<timePresolve;
    235644  handler_->printing(timeIdiot)
     
    242651  return finalStatus;
    243652}
     653// Default constructor
     654ClpSolve::ClpSolve (  )
     655{
     656  method_ = useDual;
     657  presolveType_=presolveOn;
     658  numberPasses_=5;
     659  for (int i=0;i<4;i++)
     660    options_[i]=0;
     661  for (int i=0;i<4;i++)
     662    extraInfo_[i]=-1;
     663}
     664
     665// Copy constructor.
     666ClpSolve::ClpSolve(const ClpSolve & rhs)
     667{
     668  method_ = rhs.method_;
     669  presolveType_=rhs.presolveType_;
     670  numberPasses_=rhs.numberPasses_;
     671  for (int i=0;i<4;i++)
     672    options_[i]=rhs.options_[i];
     673  for (int i=0;i<4;i++)
     674    extraInfo_[i]=rhs.extraInfo_[i];
     675}
     676// Assignment operator. This copies the data
     677ClpSolve &
     678ClpSolve::operator=(const ClpSolve & rhs)
     679{
     680  if (this != &rhs) {
     681    method_ = rhs.method_;
     682    presolveType_=rhs.presolveType_;
     683    numberPasses_=rhs.numberPasses_;
     684    for (int i=0;i<4;i++)
     685      options_[i]=rhs.options_[i];
     686    for (int i=0;i<4;i++)
     687      extraInfo_[i]=rhs.extraInfo_[i];
     688  }
     689  return *this;
     690
     691}
     692// Destructor
     693ClpSolve::~ClpSolve (  )
     694{
     695}
     696/*   which translation is:
     697     which:
     698     0 - startup in Dual  (nothing if basis exists).:
     699             0 - no basis, 1 crash
     700     1 - startup in Primal (nothing if basis exists):
     701        0 - use initiative
     702        1 - use crash
     703        2 - use idiot and look at further info
     704        3 - use sprint and look at further info
     705        4 - use all slack
     706        5 - use initiative but no idiot
     707        6 - use initiative but no sprint
     708        7 - use initiative but no crash
     709        8 - do allslack or idiot
     710        9 - do allslack or sprint
     711     2 - interrupt handling - 0 yes, 1 no (for threadsafe)
     712     3 - whether to make +- 1matrix - 0 yes, 1 no
     713*/
     714void
     715ClpSolve::setSpecialOption(int which,int value,int extraInfo)
     716{
     717  options_[which]=value;
     718  extraInfo_[which]=extraInfo;
     719}
     720int
     721ClpSolve::getSpecialOption(int which) const
     722{
     723  return options_[which];
     724}
     725
     726// Solve types
     727void
     728ClpSolve::setSolveType(SolveType method, int extraInfo)
     729{
     730  method_=method;
     731}
     732
     733ClpSolve::SolveType
     734ClpSolve::getSolveType()
     735{
     736  return method_;
     737}
     738
     739// Presolve types
     740void
     741ClpSolve::setPresolveType(PresolveType amount, int extraInfo)
     742{
     743  presolveType_=amount;
     744  numberPasses_=extraInfo;
     745}
     746ClpSolve::PresolveType
     747ClpSolve::getPresolveType()
     748{
     749  return presolveType_;
     750}
     751// Extra info for idiot (or sprint)
     752int
     753ClpSolve::getExtraInfo(int which) const
     754{
     755  return extraInfo_[which];
     756}
     757int
     758ClpSolve::getPresolvePasses() const
     759{
     760  return numberPasses_;
     761}
  • branches/pre/IdiSolve.cpp

    r208 r210  
    829829        obj[DROP-1]=best;
    830830        if (test-best<drop&&(strategy&8)==0) {
    831           if ((logLevel_&2)!=0) {
     831          if ((logLevel_&8)!=0) {
    832832          printf("Exiting as drop in %d its is %g after %d iterations\n",
    833833                 DROP*checkFrequency_,test-best,iter);
     
    10291029      maxDj = maxDj/(double) DJTEST;
    10301030      if (maxDj<djExit&&iter>50) {
    1031         printf("Exiting on low dj %g after %d iterations\n",maxDj,iter);
     1031        //printf("Exiting on low dj %g after %d iterations\n",maxDj,iter);
    10321032        break;
    10331033      }
  • branches/pre/Idiot.cpp

    r208 r210  
    99#include "ClpPresolve.hpp"
    1010#include "Idiot.hpp"
    11 #include "CoinHelperFunctions.hpp"
     11#include "CoinTime.hpp"
    1212#include "CoinMessageHandler.hpp"
    1313// Redefine stuff for Clp
     
    122122  }
    123123  sum /= (double) (nnzero+1);
    124   mu_= max(1.0e-3,sum*1.0e-5);
     124  // If mu not changed then compute
     125  if (mu_==1e-4)
     126    mu_= max(1.0e-3,sum*1.0e-5);
    125127  maxIts_=20;
    126128  maxIts2_=105;
     
    343345#endif
    344346  }
     347  int numberAway=-1;
    345348  iterationTotal = lastResult.iteration;
    346349  firstInfeas=lastResult.infeas;
     
    400403#endif
    401404    }
     405    if (iteration>50&&n==numberAway&&result.infeas<1.0e-4)
     406      break; // not much happening
     407    numberAway=n;
    402408    keepinfeas = result.infeas;
    403409    lastWeighted=result.weighted;
     
    470476            (((lambdaIteration%10)==0 && smallInfeas<keepinfeas) ||
    471477             (lambdaIteration%5)==0 && 1.5*smallInfeas<keepinfeas)) {
    472            printf(
    473    " Increasing smallInfeas from %f to %f\n",smallInfeas,1.5*smallInfeas);
     478           //printf(" Increasing smallInfeas from %f to %f\n",smallInfeas,1.5*smallInfeas);
    474479           smallInfeas *= 1.5;
    475480         }
  • branches/pre/Makefile.Clp

    r208 r210  
    66# highest level of optimization the compiler supports. If want something in
    77# between then specify the exact level you want, e.g., -O1 or -O2
    8 OptLevel := -O3
    9 #OptLevel := -g
    10 #OptLevel := -O1
     8OptLevel := -g
     9OptLevel := -O1
    1110
    1211LIBNAME := Clp
     
    5453ifeq ($(OptLevel),-g)
    5554#     CXXFLAGS += -DCLP_DEBUG
    56 #    CXXFLAGS += -DPRESOLVE_SUMMARY=1 -DDEBUG_PRESOLVE -DCHECK_CONSISTENCY=1
    5755endif
    5856ifeq ($(OptLevel),-O2)
    5957#     CXXFLAGS += -DNDEBUG
    60 #    CXXFLAGS += -DPRESOLVE_SUMMARY=1 -DDEBUG_PRESOLVE -DCHECK_CONSISTENCY=1
    6158endif
    6259
  • branches/pre/Samples/Makefile

    r200 r210  
    5959#       LDFLAGS += -lhistory -lreadline -ltermcap
    6060endif
    61 LDFLAGS += -lefence
     61#if DENSE and using given libraries
     62LDFLAGS += -llapack -lblas -lg2c
     63#LDFLAGS += -lefence
    6264###############################################################################
    6365
  • branches/pre/Samples/sprint.cpp

    r181 r210  
    125125  for (iPass=0;iPass<maxPass;iPass++) {
    126126    printf("Start of pass %d\n",iPass);
    127     printf("Bug until submodel new version\n");
     127    //printf("Bug until submodel new version\n");
    128128    CoinSort_2(sort,sort+numberSort,weight);
    129129    // Create small problem
  • branches/pre/Test/ClpMain.cpp

    r208 r210  
    1919
    2020#include "ClpFactorization.hpp"
     21#include "CoinTime.hpp"
    2122#include "ClpSimplex.hpp"
     23#include "ClpSolve.hpp"
    2224#include "ClpPackedMatrix.hpp"
    2325#include "ClpPlusMinusOneMatrix.hpp"
     
    5456  DUALTOLERANCE=1,PRIMALTOLERANCE,DUALBOUND,PRIMALWEIGHT,MAXTIME,OBJSCALE,
    5557
    56   LOGLEVEL=101,MAXFACTOR,PERTVALUE,MAXITERATION,PRESOLVEPASS,IDIOT,
     58  LOGLEVEL=101,MAXFACTOR,PERTVALUE,MAXITERATION,PRESOLVEPASS,IDIOT,SPRINT,
    5759 
    5860  DIRECTION=201,DUALPIVOT,SCALING,ERRORSALLOWED,KEEPNAMES,SPARSEFACTOR,
    59   PRIMALPIVOT,PRESOLVE,CRASH,BIASLU,PERTURBATION,
     61  PRIMALPIVOT,PRESOLVE,CRASH,BIASLU,PERTURBATION,MESSAGES,
    6062 
    6163  DIRECTORY=301,IMPORT,EXPORT,RESTORE,SAVE,DUALSIMPLEX,PRIMALSIMPLEX,
    6264  MAXIMIZE,MINIMIZE,EXIT,STDIN,UNITTEST,NETLIB_DUAL,NETLIB_PRIMAL,SOLUTION,
    63   TIGHTEN,FAKEBOUND,VERSION,PLUSMINUS,NETWORK,
     65  TIGHTEN,FAKEBOUND,HELP,PLUSMINUS,NETWORK,ALLSLACK,
    6466
    6567  INVALID=1000
    6668};
     69static void printit(const char * input)
     70{
     71  int length =strlen(input);
     72  char temp[101];
     73  int i;
     74  int n=0;
     75  for (i=0;i<length;i++) {
     76    if (input[i]=='\n') {
     77      temp[n]='\0';
     78      std::cout<<temp<<std::endl;
     79      n=0;
     80    } else if (n>=65&&input[i]==' ') {
     81      temp[n]='\0';
     82      std::cout<<temp<<std::endl;
     83      n=0;
     84    } else if (n||input[i]!=' ') {
     85      temp[n++]=input[i];
     86    }
     87  }
     88  if (n) {
     89    temp[n]='\0';
     90    std::cout<<temp<<std::endl;
     91  }
     92}
    6793/// Very simple class for setting parameters
    6894class ClpItem {
     
    79105           int lower, int upper, ClpParameterType type,bool display=true);
    80106  // Other strings will be added by insert
    81   ClpItem (std::string name, std::string help, std::string defaultValue,
    82            ClpParameterType type,bool display=true);
     107  ClpItem (std::string name, std::string help, std::string firstValue,
     108           ClpParameterType type,int defaultIndex=0,bool display=true);
    83109  // Action
    84110  ClpItem (std::string name, std::string help,
    85            ClpParameterType type,bool display=true);
     111           ClpParameterType type,int indexNumber=-1,bool display=true);
    86112  /// Copy constructor.
    87113  ClpItem(const ClpItem &);
     
    126152  inline void setCurrentOption ( int value )
    127153  { currentKeyWord_=value; };
     154  /// Sets int value
     155  inline void setIntValue ( int value )
     156  { intValue_=value; };
     157  inline int intValue () const
     158  { return intValue_; };
     159  /// Sets double value
     160  inline void setDoubleValue ( double value )
     161  { doubleValue_=value; };
     162  inline double doubleValue () const
     163  { return doubleValue_; };
     164  /// Sets string value
     165  inline void setStringValue ( std::string value )
     166  { stringValue_=value; };
     167  inline std::string stringValue () const
     168  { return stringValue_; };
    128169  /// Returns 1 if matches minimum, 2 if matches less, 0 if not matched
    129170  int matches (std::string input) const;
     
    134175  inline bool displayThis() const
    135176  { return display_;};
     177  /// Set Long help
     178  inline void setLonghelp(const std::string help)
     179  {longHelp_=help;};
     180  /// Print Long help
     181  void printLongHelp() const;
     182  /// Print action and string
     183  void printString() const;
     184  /// type for classification
     185  inline int indexNumber() const
     186  { return indexNumber_;};
    136187private:
    137188  /// gutsOfConstructor
     
    164215  std::string shortHelp_;
    165216  /// Long help
    166   std::vector<std::string> longHelp_;
     217  std::string longHelp_;
    167218  /// Action
    168219  ClpParameterType action_;
     
    171222  /// Display on ?
    172223  bool display_;
     224  /// Integer parameter - current value
     225  int intValue_;
     226  /// Double parameter - current value
     227  double doubleValue_;
     228  /// String parameter - current value
     229  std::string stringValue_;
     230  /// index number to use for display purposes
     231  int indexNumber_;
    173232  //@}
    174233};
     
    194253    action_(INVALID),
    195254    currentKeyWord_(-1),
    196     display_(false)
     255    display_(false),
     256    intValue_(-1),
     257    doubleValue_(-1.0),
     258    stringValue_(""),
     259    indexNumber_(INVALID)
    197260{
    198261}
     
    209272    action_(type),
    210273    currentKeyWord_(-1),
    211     display_(display)
     274    display_(display),
     275    intValue_(-1),
     276    doubleValue_(-1.0),
     277    stringValue_(""),
     278    indexNumber_(type)
    212279{
    213280  lowerDoubleValue_ = lower;
     
    226293    action_(type),
    227294    currentKeyWord_(-1),
    228     display_(display)
     295    display_(display),
     296    intValue_(-1),
     297    doubleValue_(-1.0),
     298    stringValue_(""),
     299    indexNumber_(type)
    229300{
    230301  gutsOfConstructor();
     
    234305// Other strings will be added by append
    235306ClpItem::ClpItem (std::string name, std::string help,
    236                   std::string defaultValue,
    237                   ClpParameterType type,bool display)
     307                  std::string firstValue,
     308                  ClpParameterType type,
     309                  int defaultIndex,bool display)
    238310  : type_(type),
    239311    lowerDoubleValue_(0.0),
     
    246318    longHelp_(),
    247319    action_(type),
    248     currentKeyWord_(0),
    249     display_(display)
     320    currentKeyWord_(defaultIndex),
     321    display_(display),
     322    intValue_(-1),
     323    doubleValue_(-1.0),
     324    stringValue_(""),
     325    indexNumber_(type)
    250326{
    251327  gutsOfConstructor();
    252   definedKeyWords_.push_back(defaultValue);
     328  definedKeyWords_.push_back(firstValue);
    253329}
    254330// Action
    255331ClpItem::ClpItem (std::string name, std::string help,
    256            ClpParameterType type,bool display)
     332           ClpParameterType type,int indexNumber,bool display)
    257333  : type_(type),
    258334    lowerDoubleValue_(0.0),
     
    266342    action_(type),
    267343    currentKeyWord_(-1),
    268     display_(display)
    269 {
     344    display_(display),
     345    intValue_(-1),
     346    doubleValue_(-1.0),
     347    stringValue_("")
     348{
     349  if (indexNumber<0)
     350    indexNumber_=type;
     351  else
     352    indexNumber_=indexNumber;
    270353  gutsOfConstructor();
    271354}
     
    290373  currentKeyWord_ = rhs.currentKeyWord_;
    291374  display_=rhs.display_;
     375  intValue_=rhs.intValue_;
     376  doubleValue_=rhs.doubleValue_;
     377  stringValue_=rhs.stringValue_;
     378  indexNumber_=rhs.indexNumber_;
    292379}
    293380
     
    320407    currentKeyWord_ = rhs.currentKeyWord_;
    321408    display_=rhs.display_;
     409    intValue_=rhs.intValue_;
     410    doubleValue_=rhs.doubleValue_;
     411    stringValue_=rhs.stringValue_;
     412    indexNumber_=rhs.indexNumber_;
    322413  }
    323414  return *this;
     
    436527  }
    437528}
     529// Print action and string
     530void
     531ClpItem::printString() const
     532{
     533  if (name_=="directory")
     534    std::cout<<"Current working directory is "<<stringValue_<<std::endl;
     535  else
     536    std::cout<<"Current default (if $ as parameter) for "<<name_
     537             <<" is "<<stringValue_<<std::endl;
     538}
    438539int
    439540ClpItem::setDoubleParameter (ClpSimplex * model,double value) const
     
    520621      model->factorization()->maximumPivots(value);
    521622      break;
    522     case DIRECTION:
    523       model->setOptimizationDirection(value);
    524       break;
    525623    case PERTVALUE:
    526624      model->setPerturbation(value);
     
    546644    value=model->factorization()->maximumPivots();
    547645    break;
    548   case DIRECTION:
    549     {
    550       double value2=model->optimizationDirection();
    551       if (value2>0.0)
    552         value=1;
    553       else if (value2<0.0)
    554         value=-1;
    555       else
    556         value=0;
    557     }
    558646    break;
    559647  case PERTVALUE:
     
    568656  }
    569657  return value;
     658}
     659// Print Long help
     660void
     661ClpItem::printLongHelp() const
     662{
     663  if (type_>=1&&type_<400) {
     664    if (type_<LOGLEVEL) {
     665      printf("Range of values is %g to %g\n",lowerDoubleValue_,upperDoubleValue_);
     666    } else if (type_<DIRECTION) {
     667      printf("Range of values is %d to %d\n",lowerIntValue_,upperIntValue_);
     668    } else if (type_<DIRECTORY) {
     669      printOptions();
     670    }
     671    printit(longHelp_.c_str());
     672  }
    570673}
    571674#ifdef READLINE     
     
    747850  {
    748851    double time1 = CoinCpuTime(),time2;
     852    // Set up all non-standard stuff
     853    //int numberModels=1;
     854    ClpSimplex * models = new ClpSimplex[1];
     855   
     856    // default action on import
     857    int allowImportErrors=0;
     858    int keepImportNames=1;
     859    int doIdiot=-2;
     860    int doCrash=0;
     861    int doSprint=-1;
     862    // set reasonable defaults
     863    int preSolve=5;
     864    models->setPerturbation(50);
     865    models->messageHandler()->setPrefix(false);
     866    std::string directory ="./";
     867    std::string importFile ="";
     868    std::string exportFile ="default.mps";
     869    std::string saveFile ="default.prob";
     870    std::string restoreFile ="default.prob";
     871    std::string solutionFile ="default.sol";
    749872#define MAXPARAMETERS 100
    750873    ClpItem parameters[MAXPARAMETERS];
    751874    int numberParameters=0;
    752875    parameters[numberParameters++]=
    753       ClpItem("?","For help",GENERALQUERY,false);
     876      ClpItem("?","For help",GENERALQUERY,-1,false);
    754877    parameters[numberParameters++]=
    755878      ClpItem("-","From stdin",
    756               STDIN,false);
     879              STDIN,299,false);
     880    parameters[numberParameters++]=
     881      ClpItem("allS!lack","Set basis back to all slack",
     882              ALLSLACK);
     883    parameters[numberParameters-1].setLonghelp
     884      (
     885       ""
     886       );
    757887    parameters[numberParameters++]=
    758888      ClpItem("biasLU","Whether factorization biased towards U",
    759               "UU",BIASLU,false);
     889              "UU",BIASLU,2,false);
    760890    parameters[numberParameters-1].append("UX");
    761891    parameters[numberParameters-1].append("LX");
     
    765895              "off",CRASH);
    766896    parameters[numberParameters-1].append("on");
     897    parameters[numberParameters-1].setLonghelp
     898      (
     899       "If crash is set on and there is an all slack basis then Clp will put structural\
     900variables into basis with the aim of getting dual feasible.  On the whole dual seems to be\
     901better without it and there alernative types of 'crash' for primal e.g. 'idiot' or 'sprint'."
     902       );
    767903    parameters[numberParameters++]=
    768904      ClpItem("direction","Minimize or Maximize",
    769905              "min!imize",DIRECTION);
    770906    parameters[numberParameters-1].append("max!imize");
    771     parameters[numberParameters++]=
    772       ClpItem("directory","Set Default import directory",
    773               DIRECTORY);
     907    parameters[numberParameters-1].append("zero");
     908    parameters[numberParameters-1].setLonghelp
     909      (
     910       "The default is minimize - use 'direction maximize' for maximization.\n\
     911You can also use the parameters 'maximize' or 'minimize'."
     912       );
     913    parameters[numberParameters++]=
     914      ClpItem("directory","Set Default directory for import etc.",
     915              DIRECTORY,299);
     916    parameters[numberParameters-1].setLonghelp
     917      (
     918       "This sets the directory which import, export, saveModel and restoreModel will use"
     919       );
     920    parameters[numberParameters-1].setStringValue(directory);
    774921    parameters[numberParameters++]=
    775922      ClpItem("dualB!ound","Initially algorithm acts as if no \
    776923gap between bounds exceeds this value",
    777924              1.0e-20,1.0e12,DUALBOUND);
     925    parameters[numberParameters-1].setLonghelp
     926      (
     927       ""
     928       );
     929    parameters[numberParameters-1].setDoubleValue(models->dualBound());
    778930    parameters[numberParameters++]=
    779931      ClpItem("dualP!ivot","Dual pivot choice algorithm",
     
    782934    parameters[numberParameters-1].append("partial");
    783935    parameters[numberParameters-1].append("steep!est");
     936    parameters[numberParameters-1].setLonghelp
     937      (
     938       ""
     939       );
    784940    parameters[numberParameters++]=
    785941      ClpItem("dualS!implex","Do dual simplex algorithm",
    786942              DUALSIMPLEX);
     943    parameters[numberParameters-1].setLonghelp
     944      (
     945       ""
     946       );
    787947    parameters[numberParameters++]=
    788948      ClpItem("dualT!olerance","For an optimal solution \
    789949no dual infeasibility may exceed this value",
    790950              1.0e-20,1.0e12,DUALTOLERANCE);
     951    parameters[numberParameters-1].setLonghelp
     952      (
     953       ""
     954       );
     955    parameters[numberParameters-1].setDoubleValue(models->dualTolerance());
    791956    parameters[numberParameters++]=
    792957      ClpItem("end","Stops clp execution",
    793958              EXIT);
     959    parameters[numberParameters-1].setLonghelp
     960      (
     961       ""
     962       );
    794963    parameters[numberParameters++]=
    795964      ClpItem("error!sAllowed","Whether to allow import errors",
    796965              "off",ERRORSALLOWED);
     966    parameters[numberParameters-1].setLonghelp
     967      (
     968       ""
     969       );
    797970    parameters[numberParameters++]=
    798971      ClpItem("exit","Stops clp execution",
    799972              EXIT);
     973    parameters[numberParameters-1].setLonghelp
     974      (
     975       ""
     976       );
    800977    parameters[numberParameters++]=
    801978      ClpItem("export","Export model as mps file",
    802979              EXPORT);
     980    parameters[numberParameters-1].setLonghelp
     981      (
     982       ""
     983       );
     984    parameters[numberParameters-1].setStringValue(exportFile);
    803985    parameters[numberParameters++]=
    804986      ClpItem("fakeB!ound","All bounds <= this value - DEBUG",
    805987              1.0,1.0e15,FAKEBOUND,false);
    806988    parameters[numberParameters++]=
     989      ClpItem("help","Print out version, non-standard options and some help",
     990              HELP);
     991    parameters[numberParameters++]=
    807992      ClpItem("idiot!Crash","Whether to try idiot crash",
    808               0,200,IDIOT);
     993              -2,200,IDIOT);
     994    parameters[numberParameters-1].setLonghelp
     995      (
     996       ""
     997       );
     998    parameters[numberParameters-1].setIntValue(doIdiot);
    809999    parameters[numberParameters++]=
    8101000      ClpItem("import","Import model from mps file",
    8111001              IMPORT);
     1002    parameters[numberParameters-1].setLonghelp
     1003      (
     1004       ""
     1005       );
     1006    parameters[numberParameters-1].setStringValue(importFile);
    8121007    parameters[numberParameters++]=
    8131008      ClpItem("keepN!ames","Whether to keep names from import",
    8141009              "on",KEEPNAMES);
    8151010    parameters[numberParameters-1].append("off");
     1011    parameters[numberParameters-1].setLonghelp
     1012      (
     1013       ""
     1014       );
    8161015    parameters[numberParameters++]=
    8171016      ClpItem("log!Level","Level of detail in output",
    8181017              0,63,LOGLEVEL);
     1018    parameters[numberParameters-1].setLonghelp
     1019      (
     1020       ""
     1021       );
     1022    parameters[numberParameters-1].setIntValue(models->logLevel());
    8191023    parameters[numberParameters++]=
    8201024      ClpItem("max!imize","Set optimization direction to maximize",
    821               MAXIMIZE);
     1025              MAXIMIZE,299);
     1026    parameters[numberParameters-1].setLonghelp
     1027      (
     1028       ""
     1029       );
    8221030    parameters[numberParameters++]=
    8231031      ClpItem("maxF!actor","Maximum number of iterations between \
    8241032refactorizations",
    8251033              1,999999,MAXFACTOR);
     1034    parameters[numberParameters-1].setLonghelp
     1035      (
     1036       ""
     1037       );
     1038    parameters[numberParameters-1].setIntValue(models->factorizationFrequency());
    8261039    parameters[numberParameters++]=
    8271040      ClpItem("maxIt!erations","Maximum number of iterations before \
    8281041stopping",
    8291042              0,99999999,MAXITERATION);
     1043    parameters[numberParameters-1].setLonghelp
     1044      (
     1045       ""
     1046       );
     1047    parameters[numberParameters-1].setIntValue(models->maximumIterations());
    8301048    parameters[numberParameters++]=
    8311049      ClpItem("min!imize","Set optimization direction to minimize",
    832               MINIMIZE);
     1050              MINIMIZE,299);
     1051    parameters[numberParameters-1].setLonghelp
     1052      (
     1053       ""
     1054       );
     1055    parameters[numberParameters++]=
     1056      ClpItem("mess!ages","Controls if Clpnnnn is printed",
     1057              "off",MESSAGES);
     1058    parameters[numberParameters-1].append("on");
     1059    parameters[numberParameters-1].setLonghelp
     1060      (""
     1061       );
    8331062    parameters[numberParameters++]=
    8341063      ClpItem("netlib","Solve entire netlib test set",
    835               NETLIB_DUAL);
     1064              NETLIB_DUAL,-1,false);
    8361065    parameters[numberParameters++]=
    8371066      ClpItem("netlibP!rimal","Solve entire netlib test set (primal)",
    838               NETLIB_PRIMAL);
     1067              NETLIB_PRIMAL,-1,false);
    8391068    parameters[numberParameters++]=
    8401069      ClpItem("network","Tries to make network matrix",
    841               NETWORK);
     1070              NETWORK,-1,false);
     1071    parameters[numberParameters-1].setLonghelp
     1072      (
     1073       ""
     1074       );
    8421075    parameters[numberParameters++]=
    8431076      ClpItem("objective!Scale","Scale factor to apply to objective",
    844               -1.0e20,1.0e20,OBJSCALE);
     1077              -1.0e20,1.0e20,OBJSCALE,false);
     1078    parameters[numberParameters-1].setDoubleValue(models->optimizationDirection());
    8451079    parameters[numberParameters++]=
    8461080      ClpItem("passP!resolve","How many passes in presolve",
    8471081              0,100,PRESOLVEPASS);
     1082    parameters[numberParameters-1].setLonghelp
     1083      (
     1084       ""
     1085       );
     1086    parameters[numberParameters-1].setIntValue(preSolve);
    8481087    parameters[numberParameters++]=
    8491088      ClpItem("pertV!alue","Method of perturbation",
    8501089              -5000,102,PERTVALUE,false);
     1090    parameters[numberParameters-1].setIntValue(models->perturbation());
    8511091    parameters[numberParameters++]=
    8521092      ClpItem("perturb!ation","Whether to perturb problem",
    8531093              "on",PERTURBATION);
    8541094    parameters[numberParameters-1].append("off");
     1095    parameters[numberParameters-1].setLonghelp
     1096      (
     1097       ""
     1098       );
    8551099    parameters[numberParameters++]=
    8561100      ClpItem("plus!Minus","Tries to make +- 1 matrix",
    857               PLUSMINUS);
     1101              PLUSMINUS,-1,false);
     1102    parameters[numberParameters-1].setLonghelp
     1103      (
     1104       ""
     1105       );
    8581106    parameters[numberParameters++]=
    8591107      ClpItem("presolve","Whether to presolve problem",
     
    8611109    parameters[numberParameters-1].append("off");
    8621110    parameters[numberParameters-1].append("more");
     1111    parameters[numberParameters-1].setLonghelp
     1112      (
     1113       ""
     1114       );
    8631115    parameters[numberParameters++]=
    8641116      ClpItem("primalP!ivot","Primal pivot choice algorithm",
    865               "steep!est",PRIMALPIVOT);
     1117              "auto!matic",PRIMALPIVOT);
    8661118    parameters[numberParameters-1].append("exa!ct");
    8671119    parameters[numberParameters-1].append("dant!zig");
     1120    parameters[numberParameters-1].append("part!ial");
     1121    parameters[numberParameters-1].append("steep!est");
     1122    //parameters[numberParameters-1].append("change");
     1123    parameters[numberParameters-1].setLonghelp
     1124      (
     1125       ""
     1126       );
    8681127    parameters[numberParameters++]=
    8691128      ClpItem("primalS!implex","Do primal simplex algorithm",
    8701129              PRIMALSIMPLEX);
     1130    parameters[numberParameters-1].setLonghelp
     1131      (
     1132       ""
     1133       );
    8711134    parameters[numberParameters++]=
    8721135      ClpItem("primalT!olerance","For an optimal solution \
    8731136no primal infeasibility may exceed this value",
    8741137              1.0e-20,1.0e12,PRIMALTOLERANCE);
     1138    parameters[numberParameters-1].setLonghelp
     1139      (
     1140       ""
     1141       );
     1142    parameters[numberParameters-1].setDoubleValue(models->primalTolerance());
    8751143    parameters[numberParameters++]=
    8761144      ClpItem("primalW!eight","Initially algorithm acts as if it \
    8771145costs this much to be infeasible",
    8781146              1.0e-20,1.0e12,PRIMALWEIGHT);
     1147    parameters[numberParameters-1].setLonghelp
     1148      (
     1149       ""
     1150       );
     1151    parameters[numberParameters-1].setDoubleValue(models->infeasibilityCost());
    8791152    parameters[numberParameters++]=
    8801153      ClpItem("quit","Stops clp execution",
    8811154              EXIT);
     1155    parameters[numberParameters-1].setLonghelp
     1156      (
     1157       ""
     1158       );
    8821159    parameters[numberParameters++]=
    8831160      ClpItem("restore!Model","Restore model from binary file",
    8841161              RESTORE);
     1162    parameters[numberParameters-1].setLonghelp
     1163      (
     1164       ""
     1165       );
     1166    parameters[numberParameters-1].setStringValue(restoreFile);
    8851167    parameters[numberParameters++]=
    8861168      ClpItem("save!Model","Save model to binary file",
    8871169              SAVE);
     1170    parameters[numberParameters-1].setLonghelp
     1171      (
     1172       ""
     1173       );
     1174    parameters[numberParameters-1].setStringValue(saveFile);
    8881175    parameters[numberParameters++]=
    8891176      ClpItem("scal!ing","Whether to scale problem",
    8901177              "on",SCALING);
    8911178    parameters[numberParameters-1].append("off");
     1179    parameters[numberParameters-1].setLonghelp
     1180      (
     1181       ""
     1182       );
    8921183    parameters[numberParameters++]=
    8931184      ClpItem("sec!onds","maximum seconds",
    8941185              0.0,1.0e12,MAXTIME);
     1186    parameters[numberParameters-1].setLonghelp
     1187      (
     1188       ""
     1189       );
     1190    parameters[numberParameters-1].setDoubleValue(models->maximumSeconds());
    8951191    parameters[numberParameters++]=
    8961192      ClpItem("sol!ution","Prints solution to file",
    8971193              SOLUTION);
     1194    parameters[numberParameters-1].setLonghelp
     1195      (
     1196       ""
     1197       );
     1198    parameters[numberParameters-1].setStringValue(solutionFile);
    8981199    parameters[numberParameters++]=
    8991200      ClpItem("spars!eFactor","Whether factorization treated as sparse",
    900               "on",SPARSEFACTOR,false);
     1201              "on",SPARSEFACTOR,0,false);
    9011202    parameters[numberParameters-1].append("off");
    9021203    parameters[numberParameters-1].append("on");
    9031204    parameters[numberParameters++]=
     1205      ClpItem("sprint!Crash","Whether to try sprint crash",
     1206              0,200,SPRINT);
     1207    parameters[numberParameters-1].setLonghelp
     1208      (
     1209       ""
     1210       );
     1211    parameters[numberParameters-1].setIntValue(doSprint);
    9041212      ClpItem("stdin","From stdin",
    905               STDIN,false);
     1213              STDIN,-1,false);
    9061214    parameters[numberParameters++]=
    9071215      ClpItem("stop","Stops clp execution",
    9081216              EXIT);
     1217    parameters[numberParameters-1].setLonghelp
     1218      (
     1219       ""
     1220       );
    9091221    parameters[numberParameters++]=
    9101222      ClpItem("tight!en","Poor person's preSolve for now",
    911               TIGHTEN);
     1223              TIGHTEN,-1,false);
    9121224    parameters[numberParameters++]=
    9131225      ClpItem("unitTest","Do unit test",
    914               UNITTEST);
    915     parameters[numberParameters++]=
    916       ClpItem("ver!sion","Print out version and non-standard options",
    917               VERSION);
     1226              UNITTEST,-1,false);
    9181227    assert(numberParameters<MAXPARAMETERS);
    9191228   
    9201229    // total number of commands read
    9211230    int numberGoodCommands=0;
    922     //int numberModels=1;
    923     ClpSimplex * models = new ClpSimplex[1];
    9241231    bool * goodModels = new bool[1];
    9251232    int getNewMatrix=0;
    9261233   
    927     // default action on import
    928     int allowImportErrors=0;
    929     int keepImportNames=1;
    930     int doIdiot=-2;
    9311234   
    9321235    int iModel=0;
    9331236    goodModels[0]=false;
    934     // set reasonable defaults
    935     int preSolve=5;
    936     models->setPerturbation(50);
    9371237    //models[0].scaling(1);
    9381238    //models[0].setDualBound(1.0e6);
     
    9431243    //ClpPrimalColumnSteepest steepP;
    9441244    //models[0].setPrimalColumnPivotAlgorithm(steepP);
    945     std::string directory ="./";
    9461245    std::string field;
    9471246    std::cout<<"Coin LP version "<<CLPVERSION
     
    9621261            <<"Clp takes input from arguments ( - switches to stdin)"
    9631262            <<std::endl
    964             <<"Enter ? for list of commands, (-)unitTest or (-)netlib"
    965             <<" for tests"<<std::endl;
     1263            <<"Enter ? for list of commands or help"<<std::endl;
    9661264          field="-";
    9671265        } else {
     
    9861284      int iParam;
    9871285      int numberMatches=0;
     1286      int firstMatch=-1;
    9881287      for ( iParam=0; iParam<numberParameters; iParam++ ) {
    9891288        int match = parameters[iParam].matches(field);
    9901289        if (match==1) {
    9911290          numberMatches = 1;
     1291          firstMatch=iParam;
    9921292          break;
    9931293        } else {
     1294          if (match&&firstMatch<0)
     1295            firstMatch=iParam;
    9941296          numberMatches += match>>1;
    9951297        }
     
    10081310          std::cout<<"abcd?? adds explanation, if only one fuller help(LATER)"<<std::endl;
    10091311          std::cout<<"abcd without value (where expected) gives current value"<<std::endl;
    1010           std::cout<<"abcd value or abcd = value sets value"<<std::endl;
     1312          std::cout<<"abcd value sets value"<<std::endl;
    10111313          std::cout<<"Commands are:"<<std::endl;
    1012           int across=0;
    1013           int maxAcross=4;
    1014           for ( iParam=0; iParam<numberParameters; iParam++ ) {
    1015             if (parameters[iParam].displayThis()) {
    1016               std::cout<<parameters[iParam].matchName()<<"  ";
    1017               across++;
    1018               if (across==maxAcross) {
    1019                 std::cout<<std::endl;
    1020                 across=0;
     1314          int maxAcross=5;
     1315          int limits[]={1,101,201,301,401};
     1316          std::vector<std::string> types;
     1317          types.push_back("Double parameters:");
     1318          types.push_back("Int parameters:");
     1319          types.push_back("Keyword parameters and others:");
     1320          types.push_back("Actions:");
     1321          int iType;
     1322          for (iType=0;iType<4;iType++) {
     1323            int across=0;
     1324            std::cout<<types[iType]<<std::endl;
     1325            for ( iParam=0; iParam<numberParameters; iParam++ ) {
     1326              int type = parameters[iParam].indexNumber();
     1327              if (parameters[iParam].displayThis()&&type>=limits[iType]
     1328                  &&type<limits[iType+1]) {
     1329                if (!across)
     1330                  std::cout<<"  ";
     1331                std::cout<<parameters[iParam].matchName()<<"  ";
     1332                across++;
     1333                if (across==maxAcross) {
     1334                  std::cout<<std::endl;
     1335                  across=0;
     1336                }
    10211337              }
    10221338            }
     1339            if (across)
     1340              std::cout<<std::endl;
    10231341          }
    1024           if (across)
    1025             std::cout<<std::endl;
    10261342        } else if (type<101) {
    10271343          // get next field as double
    10281344          double value = getDoubleField(argc,argv,&valid);
    10291345          if (!valid) {
     1346            parameters[iParam].setDoubleValue(value);
    10301347            parameters[iParam].setDoubleParameter(models+iModel,value);
    10311348          } else if (valid==1) {
     
    10331350          } else {
    10341351            std::cout<<parameters[iParam].name()<<" has value "<<
    1035               parameters[iParam].doubleParameter(models+iModel)<<std::endl;
     1352              parameters[iParam].doubleValue()<<std::endl;
    10361353          }
    10371354        } else if (type<201) {
     
    10391356          int value = getIntField(argc,argv,&valid);
    10401357          if (!valid) {
     1358            parameters[iParam].setIntValue(value);
    10411359            if (parameters[iParam].type()==PRESOLVEPASS)
    10421360              preSolve = value;
    10431361            else if (parameters[iParam].type()==IDIOT)
    10441362              doIdiot = value;
     1363            else if (parameters[iParam].type()==SPRINT)
     1364              doSprint = value;
    10451365            else
    10461366              parameters[iParam].setIntParameter(models+iModel,value);
     
    10491369          } else {
    10501370            std::cout<<parameters[iParam].name()<<" has value "<<
    1051               parameters[iParam].intParameter(models+iModel)<<std::endl;
     1371              parameters[iParam].intValue()<<std::endl;
    10521372          }
    10531373        } else if (type<301) {
     
    10711391              if (action==0)
    10721392                models[iModel].setOptimizationDirection(1);
     1393              else if (action==1)
     1394                models[iModel].setOptimizationDirection(-1);
    10731395              else
    1074                 models[iModel].setOptimizationDirection(-11);
     1396                models[iModel].setOptimizationDirection(0);
    10751397              break;
    10761398            case DUALPIVOT:
     
    10921414            case PRIMALPIVOT:
    10931415              if (action==0) {
    1094                 ClpPrimalColumnSteepest steep(1);
     1416                ClpPrimalColumnSteepest steep(3);
    10951417                models[iModel].setPrimalColumnPivotAlgorithm(steep);
    10961418              } else if (action==1) {
    10971419                ClpPrimalColumnSteepest steep(0);
    10981420                models[iModel].setPrimalColumnPivotAlgorithm(steep);
    1099               } else {
     1421              } else if (action==2) {
    11001422                ClpPrimalColumnDantzig dantzig;
    11011423                models[iModel].setPrimalColumnPivotAlgorithm(dantzig);
     1424              } else if (action==3) {
     1425                ClpPrimalColumnSteepest steep(2);
     1426                models[iModel].setPrimalColumnPivotAlgorithm(steep);
     1427              } else if (action==4) {
     1428                ClpPrimalColumnSteepest steep(0);
     1429                models[iModel].setPrimalColumnPivotAlgorithm(steep);
     1430              } else if (action==5) {
     1431                ClpPrimalColumnSteepest steep(4);
     1432                models[iModel].setPrimalColumnPivotAlgorithm(steep);
    11021433              }
    11031434              break;
     
    11321463              break;
    11331464            case CRASH:
    1134               doIdiot=-1;
     1465              doCrash=-action;
     1466              break;
     1467            case MESSAGES:
     1468              models[iModel].messageHandler()->setPrefix(action!=0);
    11351469              break;
    11361470            default:
     
    11461480          case PRIMALSIMPLEX:
    11471481            if (goodModels[iModel]) {
    1148               ClpSimplex::SolveType method;
    1149               ClpSimplex::PresolveType presolveType;
     1482              ClpSolve::SolveType method;
     1483              ClpSolve::PresolveType presolveType;
    11501484              ClpSimplex * model2 = models+iModel;
    1151               if (preSolve>5)
    1152                 presolveType=ClpSimplex::presolveMaximum;
     1485              ClpSolve solveOptions;
     1486              if (preSolve!=5&&preSolve)
     1487                presolveType=ClpSolve::presolveNumber;
    11531488              else if (preSolve)
    1154                 presolveType=ClpSimplex::presolveOn;
     1489                presolveType=ClpSolve::presolveOn;
    11551490              else
    1156                 presolveType=ClpSimplex::presolveOff;
     1491                presolveType=ClpSolve::presolveOff;
     1492              solveOptions.setPresolveType(presolveType,preSolve);
    11571493              if (type==DUALSIMPLEX)
    1158                 method=ClpSimplex::useDual;
     1494                method=ClpSolve::useDual;
    11591495              else
    1160                 method=ClpSimplex::usePrimal;
    1161               int option=0;
    1162               if (model2->perturbation()==100||method==ClpSimplex::usePrimal)
    1163                 option |= 1;
    1164               if (!model2->scalingFlag())
    1165                 option |= 2;
    1166               if (doIdiot==-1)
    1167                 option |=4;
    1168               if (doIdiot==0)
    1169                 option |= 8;
    1170               model2->initialSolve(method,presolveType,option);
     1496                method=ClpSolve::usePrimalorSprint;
     1497              solveOptions.setSolveType(method);
     1498              if (method==ClpSolve::useDual) {
     1499                // dual
     1500                if (doCrash)
     1501                  solveOptions.setSpecialOption(0,1); // crash
     1502              } else {
     1503                // primal
     1504                if (doCrash) {
     1505                  solveOptions.setSpecialOption(0,1); // crash
     1506                } else if (doSprint>0) {
     1507                  // sprint overrides idiot
     1508                  solveOptions.setSpecialOption(1,3,doSprint); // sprint
     1509                } else if (doIdiot>0) {
     1510                  solveOptions.setSpecialOption(1,2,doIdiot); // idiot
     1511                } else {
     1512                  if (doIdiot==0) {
     1513                    if (doSprint==0)
     1514                      solveOptions.setSpecialOption(1,4); // all slack
     1515                    else
     1516                      solveOptions.setSpecialOption(1,9); // all slack or sprint
     1517                  } else {
     1518                    if (doSprint==0)
     1519                      solveOptions.setSpecialOption(1,8); // all slack or idiot
     1520                    else
     1521                      solveOptions.setSpecialOption(1,7); // initiative
     1522                  }
     1523                }
     1524              }
     1525              model2->initialSolve(solveOptions);
    11711526            } else {
    11721527              std::cout<<"** Current model not valid"<<std::endl;
     
    11921547              // get next field
    11931548              field = getString(argc,argv);
     1549              if (field=="$") {
     1550                field = parameters[iParam].stringValue();
     1551              } else if (field=="EOL") {
     1552                parameters[iParam].printString();
     1553                break;
     1554              } else {
     1555                parameters[iParam].setStringValue(field);
     1556              }
    11941557              std::string fileName;
    11951558              bool canOpen=false;
     
    12351598              // get next field
    12361599              field = getString(argc,argv);
     1600              if (field=="$") {
     1601                field = parameters[iParam].stringValue();
     1602              } else if (field=="EOL") {
     1603                parameters[iParam].printString();
     1604                break;
     1605              } else {
     1606                parameters[iParam].setStringValue(field);
     1607              }
    12371608              std::string fileName;
    12381609              bool canOpen=false;
     
    12901661                    rowNames[iRow] =
    12911662                      strdup(model2->rowName(iRow).c_str());
     1663#ifdef STRIPBLANKS
     1664                    char * xx = rowNames[iRow];
     1665                    int i;
     1666                    int length = strlen(xx);
     1667                    int n=0;
     1668                    for (i=0;i<length;i++) {
     1669                      if (xx[i]!=' ')
     1670                        xx[n++]=xx[i];
     1671                    }
     1672                    xx[n]='\0';
     1673#endif
    12921674                  }
    12931675                 
     
    12961678                    columnNames[iColumn] =
    12971679                      strdup(model2->columnName(iColumn).c_str());
     1680#ifdef STRIPBLANKS
     1681                    char * xx = columnNames[iColumn];
     1682                    int i;
     1683                    int length = strlen(xx);
     1684                    int n=0;
     1685                    for (i=0;i<length;i++) {
     1686                      if (xx[i]!=' ')
     1687                        xx[n++]=xx[i];
     1688                    }
     1689                    xx[n]='\0';
     1690#endif
    12981691                  }
    12991692                }
     
    13311724              // get next field
    13321725              field = getString(argc,argv);
     1726              if (field=="$") {
     1727                field = parameters[iParam].stringValue();
     1728              } else if (field=="EOL") {
     1729                parameters[iParam].printString();
     1730                break;
     1731              } else {
     1732                parameters[iParam].setStringValue(field);
     1733              }
    13331734              std::string fileName;
    13341735              bool canOpen=false;
     
    13931794              // get next field
    13941795              field = getString(argc,argv);
     1796              if (field=="$") {
     1797                field = parameters[iParam].stringValue();
     1798              } else if (field=="EOL") {
     1799                parameters[iParam].printString();
     1800                break;
     1801              } else {
     1802                parameters[iParam].setStringValue(field);
     1803              }
    13951804              std::string fileName;
    13961805              bool canOpen=false;
     
    14271836            models[iModel].setOptimizationDirection(1);
    14281837            break;
     1838          case ALLSLACK:
     1839            models[iModel].createStatus();
     1840            break;
    14291841          case DIRECTORY:
    1430             directory = getString(argc,argv);
     1842            {
     1843              std::string name = getString(argc,argv);
     1844              if (name!="EOL") {
     1845                int length=name.length();
     1846                if (name[length-1]=='/')
     1847                  directory=name;
     1848                else
     1849                  directory = name+"/";
     1850                parameters[iParam].setStringValue(directory);
     1851              } else {
     1852                parameters[iParam].printString();
     1853              }
     1854            }
    14311855            break;
    14321856          case STDIN:
     
    15071931            }
    15081932            break;
    1509           case VERSION:
     1933          case HELP:
    15101934            std::cout<<"Coin LP version "<<CLPVERSION
    15111935                     <<", build "<<__DATE__<<std::endl;
    15121936            std::cout<<"Non default values:-"<<std::endl;
    1513             std::cout<<"Perturbation "<<models[0].perturbation()
    1514                      <<" (default 100), Presolve being done with 5 passes"
     1937            std::cout<<"Perturbation "<<models[0].perturbation()<<" (default 100)"
    15151938                     <<std::endl;
    1516             std::cout <<"Dual steepest edge steep/partial on matrix shape"
    1517                       <<std::endl;
    1518             std::cout <<"If Factorization frequency default then done on size of matrix"
    1519                       <<std::endl;
     1939            printit(
     1940                    "Presolve being done with 5 passes\n\
     1941Dual steepest edge steep/partial on matrix shape and factorization density\n\
     1942Clpnnnn taken out of messages\n\
     1943If Factorization frequency default then done on size of matrix\n\n\
     1944(-)unitTest, (-)netlib or (-)netlibp will do standard tests\n\n\
     1945You can switch to interactive mode at any time so\n\
     1946clp watson.mps -scaling off -primalsimplex\nis the same as\n\
     1947clp watson.mps -\nscaling off\nprimalsimplex"
     1948                    );
    15201949            break;
    15211950          case SOLUTION:
     
    15231952              // get next field
    15241953              field = getString(argc,argv);
     1954              if (field=="$") {
     1955                field = parameters[iParam].stringValue();
     1956              } else if (field=="EOL") {
     1957                parameters[iParam].printString();
     1958                break;
     1959              } else {
     1960                parameters[iParam].setStringValue(field);
     1961              }
    15251962              std::string fileName;
    15261963              FILE *fp=NULL;
     
    16032040      } else if (numberMatches==1) {
    16042041        if (!numberQuery) {
    1605           std::cout<<"Short match for "<<field<<" possible completion:"
    1606                    <<std::endl;
    1607           for ( iParam=0; iParam<numberParameters; iParam++ ) {
    1608             int match = parameters[iParam].matches(field);
    1609             if (match&&parameters[iParam].displayThis())
    1610               std::cout<<parameters[iParam].matchName()<<std::endl;
    1611           }
     2042          std::cout<<"Short match for "<<field<<" - completion: ";
     2043          std::cout<<parameters[firstMatch].matchName()<<std::endl;
    16122044        } else if (numberQuery) {
    1613           std::cout<<"Short match for "<<field<<" completion:"
    1614                    <<std::endl;
    1615           for ( iParam=0; iParam<numberParameters; iParam++ ) {
    1616             int match = parameters[iParam].matches(field);
    1617             if (match&&parameters[iParam].displayThis()) {
    1618               std::cout<<parameters[iParam].matchName()<<" : ";
    1619               std::cout<<parameters[iParam].shortHelp()<<std::endl;
    1620             }
    1621           }
     2045          std::cout<<parameters[firstMatch].matchName()<<" : ";
     2046          std::cout<<parameters[firstMatch].shortHelp()<<std::endl;
     2047          if (numberQuery>=2)
     2048            parameters[firstMatch].printLongHelp();
    16222049        }
    16232050      } else {
  • branches/pre/Test/unitTest.cpp

    r208 r210  
    1414#include "CoinPackedVector.hpp"
    1515#include "CoinHelperFunctions.hpp"
     16#include "CoinTime.hpp"
    1617
    1718#include "ClpFactorization.hpp"
  • branches/pre/include/ClpMatrixBase.hpp

    r205 r210  
    5858    /** Delete the rows whose indices are listed in <code>indDel</code>. */
    5959    virtual void deleteRows(const int numDel, const int * indDel) = 0;
     60  /// Append Columns
     61  virtual void appendCols(int number, const CoinPackedVectorBase * const * columns);
     62  /// Append Rows
     63  virtual void appendRows(int number, const CoinPackedVectorBase * const * rows);
    6064
    6165  /** Returns a new matrix in reverse order without gaps
  • branches/pre/include/ClpMessage.hpp

    r208 r210  
    6868  CLP_TIMING,
    6969  CLP_INTERVAL_TIMING,
     70  CLP_SPRINT,
    7071  CLP_DUMMY_END
    7172};
  • branches/pre/include/ClpModel.hpp

    r205 r210  
    220220  inline void setProblemStatus(int problemStatus)
    221221  { problemStatus_ = problemStatus;};
     222   /** Secondary status of problem - may get extended
     223       0 - none
     224       1 - primal infeasible because dual limit reached
     225   */
     226   inline int secondaryStatus() const            { return secondaryStatus_; }
     227  inline void setSecondaryStatus(int status)
     228  { secondaryStatus_ = status;};
    222229   /// Are there a numerical difficulties?
    223230   bool isAbandoned() const             { return problemStatus_==4; }
     
    521528  /// Status of problem
    522529  int problemStatus_;
     530  /// Secondary status of problem
     531  int secondaryStatus_;
    523532  /// length of names (0 means no names)
    524533  int lengthNames_;
  • branches/pre/include/ClpNonLinearCost.hpp

    r185 r210  
    100100   May need to be inline for speed */
    101101  double setOne(int sequence, double solutionValue);
     102  /** Sets bounds and cost for outgoing variable
     103      may change value
     104      Returns direction */
     105  int setOneOutgoing(int sequence, double &solutionValue);
    102106  /// Returns nearest bound
    103107  double nearest(int sequence, double solutionValue);
  • branches/pre/include/ClpPackedMatrix.hpp

    r205 r210  
    5858    virtual void deleteRows(const int numDel, const int * indDel)
    5959  { matrix_->deleteRows(numDel,indDel);};
     60  /// Append Columns
     61  virtual void appendCols(int number, const CoinPackedVectorBase * const * columns)
     62  { matrix_->appendCols(number,columns);};
     63  /// Append Rows
     64  virtual void appendRows(int number, const CoinPackedVectorBase * const * rows)
     65  { matrix_->appendRows(number,rows);};
    6066  /** Replace the elements of a vector.  The indices remain the same.
    6167      This is only needed if scaling and a row copy is used.
  • branches/pre/include/ClpPlusMinusOneMatrix.hpp

    r207 r210  
    5050    /** Delete the rows whose indices are listed in <code>indDel</code>. */
    5151  virtual void deleteRows(const int numDel, const int * indDel);
     52  /// Append Columns
     53  virtual void appendCols(int number, const CoinPackedVectorBase * const * columns);
     54  /// Append Rows
     55  virtual void appendRows(int number, const CoinPackedVectorBase * const * rows);
    5256  /** Returns a new matrix in reverse order without gaps */
    5357  virtual ClpMatrixBase * reverseOrderedCopy() const;
  • branches/pre/include/ClpPrimalColumnSteepest.hpp

    r180 r210  
    6565  ///@name Constructors and destructors
    6666  //@{
    67   /// Default Constructor
    68   ClpPrimalColumnSteepest(int mode=0);
     67  /** Default Constructor
     68      0 is exact devex, 1 full steepest, 2 is partial exact devex
     69      3 switches between 0 and 2 depending on factorization
     70      4 starts as partial dantzig but then may switch between 0 and 2.
     71      By partial exact devex is meant that the weights are updated as normal
     72      but only part of the nonbasic variables are scanned. 
     73      This can be faster on very easy problems.
     74  */
     75  ClpPrimalColumnSteepest(int mode=3);
    6976 
    7077  /// Copy constructor
     
    120127  */
    121128  int state_;
    122   /// If 0 then we are using exact devex, 1 then full
     129  /**
     130      0 is exact devex, 1 full steepest, 2 is partial exact devex
     131      3 switches between 0 and 2 depending on factorization
     132      4 starts as partial dantzig but then may switch between 0 and 2.
     133      By partial exact devex is meant that the weights are updated as normal
     134      but only part of the nonbasic variables are scanned. 
     135      This can be faster on very easy problems.*/
    123136  int mode_;
     137  /// Number of times switched from partial dantzig to 0/2
     138  int numberSwitched_;
    124139  // This is pivot row (or pivot sequence round re-factorization)
    125140  int pivotSequence_; 
  • branches/pre/include/ClpSimplex.hpp

    r208 r210  
    1515#include "ClpModel.hpp"
    1616#include "ClpMatrixBase.hpp"
     17#include "ClpSolve.hpp"
    1718class ClpDualRowPivot;
    1819class ClpPrimalColumnPivot;
     
    4849
    4950public:
    50 
    51   /** enums for solve function */
    52   enum SolveType {
    53     useDual=0,
    54     usePrimal,
    55     useSprint,
    56     automatic
    57   };
    58   enum PresolveType {
    59     presolveOn=0,
    60     presolveOff,
    61     presolveMaximum
    62   };
    6351  /** enums for status of various sorts.
    6452      First 4 match CoinWarmStartBasis,
     
    155143  /**@name Functions most useful to user */
    156144  //@{
    157   /** General solve algorithm which can do presolve
    158       special options (bits)
    159       1 - do not perturb
    160       2 - do not scale
    161       4 - use crash (default allslack in dual, idiot in primal)
    162       8 - all slack basis in primal
    163       16 - switch off interrupt handling
    164       32 - do not try and make plus minus one matrix
     145  /** General solve algorithm which can do presolve.
     146      See  ClpSolve.hpp for options
    165147   */
    166   int initialSolve(SolveType method=useDual, PresolveType presolve=presolveOn,
    167                    int specialOptions=0);
     148  int initialSolve(ClpSolve & options);
    168149  /** Dual algorithm - see ClpSimplexDual.hpp for method */
    169150  int dual(int ifValuesPass=0);
     
    601582    st_byte |= status;
    602583  };
     584  /** Normally the first factorization does sparse coding because
     585      the factorization could be singular.  This allows initial dense
     586      factorization when it is known to be safe
     587  */
     588  void setInitialDenseFactorization(bool onOff);
     589  bool  initialDenseFactorization() const;
    603590  /** Return sequence In or Out */
    604591  inline int sequenceIn() const
     
    909896      2 - Keep nonLinearCost round solves
    910897      4 - Force outgoing variables to exact bound (primal)
     898      8 - Safe to use dense initial factorization
    911899  */
    912900  unsigned int specialOptions_;
  • branches/pre/include/ClpSimplexPrimal.hpp

    r180 r210  
    212212                             ClpSimplexProgress * progress);
    213213  /// Perturbs problem (method depends on perturbation())
    214   void perturb();
    215   /// Take off effect of perturbation
    216   void unPerturb();
     214  void perturb(int type);
     215  /// Take off effect of perturbation and say whether to try dual
     216  bool unPerturb();
    217217  /// Unflag all variables and return number unflagged
    218218  int unflag();
Note: See TracChangeset for help on using the changeset viewer.