Changeset 1111


Ignore:
Timestamp:
Sep 21, 2007 4:13:41 PM (12 years ago)
Author:
forrest
Message:

trying to improve networks

Location:
trunk/Clp
Files:
1 added
14 edited

Legend:

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

    r1108 r1111  
    14941494  parameters[numberParameters++]=
    14951495    CbcOrClpParam("dualize","Solves dual reformulation",
    1496                   0,2,DUALIZE,false);
     1496                  0,3,DUALIZE,false);
    14971497  parameters[numberParameters-1].setLonghelp
    14981498    (
  • trunk/Clp/src/ClpFactorization.cpp

    r919 r1111  
    278278      if (networkMatrix)
    279279        biasLU_=0; // All to U if network
    280       int saveMaximumPivots = maximumPivots();
     280      //int saveMaximumPivots = maximumPivots();
    281281      delete networkBasis_;
    282282      networkBasis_ = NULL;
     
    414414#ifndef SLIM_CLP
    415415        if (networkMatrix) {
    416           maximumPivots(saveMaximumPivots);
     416          maximumPivots(CoinMax(2000,maximumPivots()));
     417          // redo arrays
     418          for (int iRow=0;iRow<4;iRow++) {
     419            int length =model->numberRows()+maximumPivots();
     420            if (iRow==3||model->objectiveAsObject()->type()>1)
     421              length += model->numberColumns();
     422            model->rowArray(iRow)->reserve(length);
     423          }
    417424          // create network factorization
    418425          if (doCheck)
     
    614621    // network - fake factorization - do nothing
    615622    status_=0;
     623    numberPivots_ = 0;
    616624  }
    617625#endif
  • trunk/Clp/src/ClpMain.cpp

    r1092 r1111  
    1616#include "CoinSort.hpp"
    1717// History since 1.0 at end
    18 #define CLPVERSION "1.05.00"
     18#define CLPVERSION "1.06.00"
    1919
    2020#include "CoinMpsIO.hpp"
     
    3636#include "ClpPresolve.hpp"
    3737#include "CbcOrClpParam.hpp"
     38#include "CoinSignal.hpp"
    3839#ifdef DMALLOC
    3940#include "dmalloc.h"
     
    5253static bool maskMatches(const int * starts, char ** masks,
    5354                        std::string & check);
     55static ClpSimplex * currentModel = NULL;
     56
     57extern "C" {
     58   static void signal_handler(int whichSignal)
     59   {
     60      if (currentModel!=NULL)
     61         currentModel->setMaximumIterations(0); // stop at next iterations
     62      return;
     63   }
     64}
    5465
    5566//#############################################################################
     
    115126    int basisHasValues=0;
    116127    int substitution=3;
    117     int dualize=0;
     128    int dualize=3;  // dualize if looks promising
    118129    std::string exportBasisFile ="default.bas";
    119130    std::string saveFile ="default.prob";
     
    439450                models[iModel].setPrimalColumnPivotAlgorithm(dantzig);
    440451              } else if (action==3) {
    441                 ClpPrimalColumnSteepest steep(2);
     452                ClpPrimalColumnSteepest steep(4);
    442453                models[iModel].setPrimalColumnPivotAlgorithm(steep);
    443454              } else if (action==4) {
     
    445456                models[iModel].setPrimalColumnPivotAlgorithm(steep);
    446457              } else if (action==5) {
    447                 ClpPrimalColumnSteepest steep(4);
     458                ClpPrimalColumnSteepest steep(2);
    448459                models[iModel].setPrimalColumnPivotAlgorithm(steep);
    449460              } else if (action==6) {
     
    559570              ClpSimplex * model2 = models+iModel;
    560571              if (dualize) {
    561                 model2 = ((ClpSimplexOther *) model2)->dualOfModel();
    562                 printf("Dual of model has %d rows and %d columns\n",
    563                        model2->numberRows(),model2->numberColumns());
    564                 model2->setOptimizationDirection(1.0);
     572                bool tryIt=true;
     573                double fractionColumn=1.0;
     574                double fractionRow=1.0;
     575                if (dualize==3) {
     576                  dualize=1;
     577                  int numberColumns=model2->numberColumns();
     578                  int numberRows=model2->numberRows();
     579                  if (numberRows<50000||5*numberColumns>numberRows) {
     580                    tryIt=false;
     581                  } else {
     582                    fractionColumn=0.1;
     583                    fractionRow=0.1;
     584                  }
     585                }
     586                if (tryIt) {
     587                  model2 = ((ClpSimplexOther *) model2)->dualOfModel(fractionRow,fractionColumn);
     588                  if (model2) {
     589                    printf("Dual of model has %d rows and %d columns\n",
     590                           model2->numberRows(),model2->numberColumns());
     591                    model2->setOptimizationDirection(1.0);
     592                  } else {
     593                    model2 = models+iModel;
     594                    dualize=0;
     595                  }
     596                } else {
     597                  dualize=0;
     598                }
    565599              }
    566600              ClpSolve solveOptions;
     
    694728              if (dualize) {
    695729                int returnCode=((ClpSimplexOther *) models+iModel)->restoreFromDual(model2);
     730                if (model2->status()==3)
     731                  returnCode=0;
    696732                delete model2;
    697                 if (returnCode&&dualize!=2)
     733                if (returnCode&&dualize!=2) {
     734                  CoinSighandler_t saveSignal=SIG_DFL;
     735                  currentModel = models+iModel;
     736                  // register signal handler
     737                  saveSignal = signal(SIGINT,signal_handler);
    698738                  models[iModel].primal(1);
     739                  currentModel=NULL;
     740                }
    699741              }
    700742              if (status>=0)
     
    9781020                bool deleteModel2=false;
    9791021                ClpSimplex * model2 = models+iModel;
    980                 if (dualize) {
     1022                if (dualize&&dualize<3) {
    9811023                  model2 = ((ClpSimplexOther *) model2)->dualOfModel();
    9821024                  printf("Dual of model has %d rows and %d columns\n",
  • trunk/Clp/src/ClpMatrixBase.hpp

    r1055 r1111  
    335335                                       int numberRows, const int * whichRows,
    336336                                       int numberColumns, const int * whichColumns) const;
    337  
     337  /// Gets rid of any mutable by products
     338  virtual void backToBasics() {}
    338339  /** Returns type.
    339340      The types which code may need to know about are:
  • trunk/Clp/src/ClpNetworkMatrix.cpp

    r1034 r1111  
    88#include "CoinIndexedVector.hpp"
    99#include "CoinHelperFunctions.hpp"
     10#include "CoinPackedVector.hpp"
    1011
    1112#include "ClpSimplex.hpp"
     
    1617#include "ClpMessage.hpp"
    1718#include <iostream>
     19#include <cassert>
    1820
    1921//#############################################################################
     
    2830{
    2931  setType(11);
    30   elements_ = NULL;
    31   starts_ = NULL;
     32  matrix_ = NULL;
    3233  lengths_=NULL;
    3334  indices_=NULL;
     
    4344{
    4445  setType(11);
    45   elements_ = NULL;
    46   starts_ = NULL;
     46  matrix_ = NULL;
    4747  lengths_=NULL;
    4848  indices_=new int[2*numberColumns];;
     
    6868: ClpMatrixBase(rhs)
    6969
    70   elements_ = NULL;
    71   starts_ = NULL;
     70  matrix_ = NULL;
    7271  lengths_=NULL;
    7372  indices_=NULL;
     
    9190
    9291  setType(11);
    93   elements_ = NULL;
    94   starts_ = NULL;
     92  matrix_ = NULL;
    9593  lengths_=NULL;
    9694  indices_=NULL;
     
    187185ClpNetworkMatrix::~ClpNetworkMatrix ()
    188186{
    189   delete [] elements_;
    190   delete [] starts_;
     187  delete matrix_;
    191188  delete [] lengths_;
    192189  delete [] indices_;
     
    201198  if (this != &rhs) {
    202199    ClpMatrixBase::operator=(rhs);
    203     delete [] elements_;
    204     delete [] starts_;
     200    delete matrix_;
    205201    delete [] lengths_;
    206202    delete [] indices_;
    207     elements_ = NULL;
    208     starts_ = NULL;
     203    matrix_ = NULL;
    209204    lengths_=NULL;
    210205    indices_=NULL;
     
    669664ClpNetworkMatrix::getPackedMatrix() const
    670665{
    671   return new CoinPackedMatrix(true,numberRows_,numberColumns_,
    672                               2*numberColumns_,
    673                               getElements(),indices_,
    674                               getVectorStarts(),getVectorLengths());
    675 
     666  if (!matrix_) {
     667    assert (trueNetwork_); // fix later
     668    int numberElements = 2*numberColumns_;
     669    double * elements = new double [numberElements];
     670    CoinBigIndex i;
     671    for (i=0;i<2*numberColumns_;i+=2) {
     672      elements[i]=-1.0;
     673      elements[i+1]=1.0;
     674    }
     675    CoinBigIndex * starts = new CoinBigIndex [numberColumns_+1];
     676    for (i=0;i<numberColumns_+1;i++) {
     677      starts[i]=2*i;
     678    }
     679    // use assignMatrix to save space
     680    delete [] lengths_;
     681    lengths_ = NULL;
     682    matrix_ = new CoinPackedMatrix();
     683    int * indices = CoinCopyOfArray(indices_,2*numberColumns_);
     684    matrix_->assignMatrix(true,numberRows_,numberColumns_,
     685                          getNumElements(),
     686                          elements,indices,
     687                          starts,lengths_);
     688    assert(!elements);
     689    assert(!starts);
     690    assert (!indices);
     691    assert (!lengths_);
     692  }
     693  return matrix_;
    676694}
    677695/* A vector containing the elements in the packed matrix. Note that there
     
    682700ClpNetworkMatrix::getElements() const
    683701{
    684   assert (trueNetwork_); // fix later
    685   if (!elements_) {
    686     elements_ = new double [2*numberColumns_];
    687     int i;
    688     for (i=0;i<2*numberColumns_;i+=2) {
    689       elements_[i]=-1.0;
    690       elements_[i+1]=1.0;
    691     }
    692   }
    693   return elements_;
     702  if (!matrix_)
     703    getPackedMatrix();
     704  return matrix_->getElements();
    694705}
    695706
     
    697708ClpNetworkMatrix::getVectorStarts() const
    698709{
    699   assert (trueNetwork_); // fix later
    700   if (!starts_) {
    701     starts_ = new CoinBigIndex [numberColumns_+1];
    702     int i;
    703     for (i=0;i<numberColumns_+1;i++) {
    704       starts_[i]=2*i;
    705     }
    706   }
    707   return starts_;
     710  if (!matrix_)
     711    getPackedMatrix();
     712  return matrix_->getVectorStarts();
    708713}
    709714/* The lengths of the major-dimension vectors. */
     
    724729void ClpNetworkMatrix::deleteCols(const int numDel, const int * indDel)
    725730{
    726   std::cerr<<"deleteCols not implemented in ClpNetworkMatrix"<<std::endl;
    727   abort();
     731  assert (trueNetwork_);
     732  int iColumn;
     733  int numberBad=0;
     734  // Use array to make sure we can have duplicates
     735  char * which = new char[numberColumns_];
     736  memset(which,0,numberColumns_);
     737  int nDuplicate=0;
     738  for (iColumn=0;iColumn<numDel;iColumn++) {
     739    int jColumn = indDel[iColumn];
     740    if (jColumn<0||jColumn>=numberColumns_) {
     741      numberBad++;
     742    } else {
     743      if (which[jColumn])
     744        nDuplicate++;
     745      else
     746        which[jColumn]=1;
     747    }
     748  }
     749  if (numberBad)
     750    throw CoinError("Indices out of range", "deleteCols", "ClpNetworkMatrix");
     751  int newNumber = numberColumns_-numDel+nDuplicate;
     752  // Get rid of temporary arrays
     753  delete [] lengths_;
     754  lengths_=NULL;
     755  delete matrix_;
     756  matrix_= NULL;
     757  int newSize=2*newNumber;
     758  int * newIndices = new int [newSize];
     759  newSize=0;
     760  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
     761    if (!which[iColumn]) {
     762      CoinBigIndex start;
     763      CoinBigIndex i;
     764      start = 2*iColumn;
     765      for (i=start;i<start+2;i++)
     766        newIndices[newSize++]=indices_[i];
     767    }
     768  }
     769  delete [] which;
     770  delete [] indices_;
     771  indices_= newIndices;
     772  numberColumns_ = newNumber;
    728773}
    729774/* Delete the rows whose indices are listed in <code>indDel</code>. */
    730775void ClpNetworkMatrix::deleteRows(const int numDel, const int * indDel)
    731776{
    732   std::cerr<<"deleteRows not implemented in ClpNetworkMatrix"<<std::endl;
    733   abort();
     777  int iRow;
     778  int numberBad=0;
     779  // Use array to make sure we can have duplicates
     780  int * which = new int [numberRows_];
     781  memset(which,0,numberRows_*sizeof(int));
     782  for (iRow=0;iRow<numDel;iRow++) {
     783    int jRow = indDel[iRow];
     784    if (jRow<0||jRow>=numberRows_) {
     785      numberBad++;
     786    } else {
     787      which[jRow]=1;
     788    }
     789  }
     790  if (numberBad)
     791    throw CoinError("Indices out of range", "deleteRows", "ClpNetworkMatrix");
     792  // Only valid of all columns have 0 entries
     793  int iColumn;
     794  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
     795    CoinBigIndex start;
     796    CoinBigIndex i;
     797    start = 2*iColumn;
     798    for (i=start;i<start+2;i++) {
     799      int iRow = indices_[i];
     800      if (which[iRow])
     801        numberBad++;
     802    }
     803  }
     804  if (numberBad)
     805    throw CoinError("Row has entries", "deleteRows", "ClpNetworkMatrix");
     806  int newNumber=0;
     807  for (iRow=0;iRow<numberRows_;iRow++) {
     808    if (!which[iRow])
     809      which[iRow]=newNumber++;
     810    else
     811      which[iRow]=-1;
     812  }
     813  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
     814    CoinBigIndex start;
     815    CoinBigIndex i;
     816    start = 2*iColumn;
     817    for (i=start;i<start+2;i++) {
     818      int iRow = indices_[i];
     819      indices_[i]=which[iRow];
     820    }
     821  }
     822  delete [] which;
     823  numberRows_ = newNumber;
    734824}
    735825/* Given positive integer weights for each row fills in sum of weights
     
    10181108ClpNetworkMatrix::releasePackedMatrix() const
    10191109{
    1020   delete [] elements_;
     1110  delete matrix_;
    10211111  delete [] lengths_;
    1022   elements_=NULL;
     1112  matrix_=NULL;
    10231113  lengths_=NULL;
    10241114}
     1115// Append Columns
     1116void
     1117ClpNetworkMatrix::appendCols(int number, const CoinPackedVectorBase * const * columns)
     1118{
     1119  int iColumn;
     1120  int numberBad=0;
     1121  for (iColumn=0;iColumn<number;iColumn++) {
     1122    int n=columns[iColumn]->getNumElements();
     1123    const double * element = columns[iColumn]->getElements();
     1124    if (n!=2)
     1125      numberBad++;
     1126    if (fabs(element[0])!=1.0||fabs(element[1])!=1.0)
     1127        numberBad++;
     1128    else if (element[0]*element[1]!=-1.0)
     1129        numberBad++;
     1130  }
     1131  if (numberBad)
     1132    throw CoinError("Not network", "appendCols", "ClpNetworkMatrix");
     1133  // Get rid of temporary arrays
     1134  delete [] lengths_;
     1135  lengths_=NULL;
     1136  delete matrix_;
     1137  matrix_= NULL;
     1138  CoinBigIndex size = 2*number;
     1139  int * temp2 = new int [numberColumns_*2+size];
     1140  memcpy(temp2,indices_,numberColumns_*2*sizeof(int));
     1141  delete [] indices_;
     1142  indices_= temp2;
     1143  // now add
     1144  size=2*numberColumns_;
     1145  for (iColumn=0;iColumn<number;iColumn++) {
     1146    const int * row = columns[iColumn]->getIndices();
     1147    const double * element = columns[iColumn]->getElements();
     1148    if (element[0]==-1.0) {
     1149      indices_[size++]=row[0];
     1150      indices_[size++]=row[1];
     1151    } else {
     1152      indices_[size++]=row[1];
     1153      indices_[size++]=row[0];
     1154    }
     1155  }
     1156 
     1157  numberColumns_ += number;
     1158}
     1159// Append Rows
     1160void
     1161ClpNetworkMatrix::appendRows(int number, const CoinPackedVectorBase * const * rows)
     1162{
     1163  // must be zero arrays
     1164  int numberBad=0;
     1165  int iRow;
     1166  for (iRow=0;iRow<number;iRow++) {
     1167    numberBad += rows[iRow]->getNumElements();
     1168  }
     1169  if (numberBad)
     1170    throw CoinError("Not NULL rows", "appendRows", "ClpNetworkMatrix");
     1171  numberRows_ += number;
     1172}
     1173#ifndef SLIM_CLP
     1174/* Append a set of rows/columns to the end of the matrix. Returns number of errors
     1175   i.e. if any of the new rows/columns contain an index that's larger than the
     1176   number of columns-1/rows-1 (if numberOther>0) or duplicates
     1177   If 0 then rows, 1 if columns */
     1178int
     1179ClpNetworkMatrix::appendMatrix(int number, int type,
     1180                                    const CoinBigIndex * starts, const int * index,
     1181                                    const double * element, int numberOther)
     1182{
     1183  int numberErrors=0;
     1184  // make into CoinPackedVector
     1185  CoinPackedVectorBase ** vectors=
     1186    new CoinPackedVectorBase * [number];
     1187  int iVector;
     1188  for (iVector=0;iVector<number;iVector++) {
     1189    int iStart = starts[iVector];
     1190    vectors[iVector] =
     1191      new CoinPackedVector(starts[iVector+1]-iStart,
     1192                           index+iStart,element+iStart);
     1193  }
     1194  if (type==0) {
     1195    // rows
     1196    appendRows(number,vectors);
     1197  } else {
     1198    // columns
     1199    appendCols(number,vectors);
     1200  }
     1201  for (iVector=0;iVector<number;iVector++)
     1202    delete vectors[iVector];
     1203  delete [] vectors;
     1204  return numberErrors;
     1205}
     1206#endif
     1207/* Subset clone (without gaps).  Duplicates are allowed
     1208   and order is as given */
     1209ClpMatrixBase *
     1210ClpNetworkMatrix::subsetClone (int numberRows, const int * whichRows,
     1211                              int numberColumns,
     1212                              const int * whichColumns) const
     1213{
     1214  return new ClpNetworkMatrix(*this, numberRows, whichRows,
     1215                                   numberColumns, whichColumns);
     1216}
     1217/* Subset constructor (without gaps).  Duplicates are allowed
     1218   and order is as given */
     1219ClpNetworkMatrix::ClpNetworkMatrix (
     1220                       const ClpNetworkMatrix & rhs,
     1221                       int numberRows, const int * whichRow,
     1222                       int numberColumns, const int * whichColumn)
     1223: ClpMatrixBase(rhs)
     1224
     1225  setType(11);
     1226  matrix_ = NULL;
     1227  lengths_=NULL;
     1228  indices_=new int[2*numberColumns];;
     1229  numberRows_=numberRows;
     1230  numberColumns_=numberColumns;
     1231  trueNetwork_=true;
     1232  int iColumn;
     1233  int numberBad=0;
     1234  int * which = new int [rhs.numberRows_];
     1235  int iRow;
     1236  for (iRow=0;iRow<rhs.numberRows_;iRow++)
     1237    which[iRow]=-1;
     1238  int n=0;
     1239  for (iRow=0;iRow<numberRows;iRow++) {
     1240    int jRow=whichRow[iRow];
     1241    assert (jRow>=0&&jRow<rhs.numberRows_);
     1242    which[jRow]=n++;
     1243  }
     1244  for (iColumn=0;iColumn<numberColumns;iColumn++) {
     1245    CoinBigIndex start;
     1246    CoinBigIndex i;
     1247    start = 2*iColumn;
     1248    CoinBigIndex offset = 2*whichColumn[iColumn]-start;
     1249    for (i=start;i<start+2;i++) {
     1250      int iRow = rhs.indices_[i+offset];
     1251      iRow = which[iRow];
     1252      if (iRow<0)
     1253        numberBad++;
     1254      else
     1255        indices_[i]=iRow;
     1256    }
     1257  }
     1258  if (numberBad)
     1259    throw CoinError("Invalid rows", "subsetConstructor", "ClpNetworkMatrix");
     1260}
  • trunk/Clp/src/ClpNetworkMatrix.hpp

    r1055 r1111  
    5252    /** Delete the rows whose indices are listed in <code>indDel</code>. */
    5353  virtual void deleteRows(const int numDel, const int * indDel);
     54  /// Append Columns
     55  virtual void appendCols(int number, const CoinPackedVectorBase * const * columns);
     56  /// Append Rows
     57  virtual void appendRows(int number, const CoinPackedVectorBase * const * rows);
     58#ifndef SLIM_CLP
     59  /** Append a set of rows/columns to the end of the matrix. Returns number of errors
     60      i.e. if any of the new rows/columns contain an index that's larger than the
     61      number of columns-1/rows-1 (if numberOther>0) or duplicates
     62      If 0 then rows, 1 if columns */
     63  virtual int appendMatrix(int number, int type,
     64                           const CoinBigIndex * starts, const int * index,
     65                           const double * element, int numberOther=-1);
     66#endif
    5467  /** Returns a new matrix in reverse order without gaps */
    5568  virtual ClpMatrixBase * reverseOrderedCopy() const;
     
    171184  /// Clone
    172185  virtual ClpMatrixBase * clone() const ;
     186  /** Subset constructor (without gaps).  Duplicates are allowed
     187      and order is as given */
     188  ClpNetworkMatrix (const ClpNetworkMatrix & wholeModel,
     189                    int numberRows, const int * whichRows,
     190                    int numberColumns, const int * whichColumns);
     191  /** Subset clone (without gaps).  Duplicates are allowed
     192      and order is as given */
     193  virtual ClpMatrixBase * subsetClone (
     194                    int numberRows, const int * whichRows,
     195                    int numberColumns, const int * whichColumns) const ;
    173196   //@}
    174197   
     
    179202   //@{
    180203  /// For fake CoinPackedMatrix
    181   mutable double * elements_;
    182   mutable CoinBigIndex * starts_;
     204  mutable CoinPackedMatrix * matrix_;
    183205  mutable int * lengths_;
    184206  /// Data -1, then +1 rows in pairs (row==-1 if one entry)
  • trunk/Clp/src/ClpPlusMinusOneMatrix.cpp

    r1106 r1111  
    10621062                              startPositive_,getVectorLengths());
    10631063    delete [] elements;
     1064    delete [] lengths_;
     1065    lengths_ = NULL;
    10641066  }
    10651067  return matrix_;
     
    11801182  }
    11811183  if (numberBad)
    1182     throw CoinError("Indices out of range", "deleteCols", "ClpPlusMinusOneMatrix");
     1184    throw CoinError("Indices out of range", "deleteRows", "ClpPlusMinusOneMatrix");
    11831185  CoinBigIndex iElement;
    11841186  CoinBigIndex numberElements=startPositive_[numberColumns_];
  • trunk/Clp/src/ClpSimplex.cpp

    r1103 r1111  
    4545  columnDualInfeasibility_(0.0),
    4646  rowDualInfeasibility_(0.0),
    47   moreSpecialOptions_(0),
     47  moreSpecialOptions_(2),
    4848  rowDualSequence_(-2),
    4949  primalToleranceToGetOptimal_(-1.0),
     
    156156  columnDualInfeasibility_(0.0),
    157157  rowDualInfeasibility_(0.0),
    158   moreSpecialOptions_(0),
     158  moreSpecialOptions_(2),
    159159  rowDualSequence_(-2),
    160160  primalToleranceToGetOptimal_(-1.0),
     
    16151615  columnDualInfeasibility_(0.0),
    16161616  rowDualInfeasibility_(0.0),
    1617   moreSpecialOptions_(0),
     1617  moreSpecialOptions_(2),
    16181618  rowDualSequence_(-2),
    16191619  primalToleranceToGetOptimal_(-1.0),
     
    17191719  columnDualInfeasibility_(0.0),
    17201720  rowDualInfeasibility_(0.0),
    1721   moreSpecialOptions_(0),
     1721  moreSpecialOptions_(2),
    17221722  rowDualSequence_(-2),
    17231723  primalToleranceToGetOptimal_(-1.0),
     
    40914091  CoinPackedMatrix copy;
    40924092  copy.reverseOrderedCopyOf(*matrix());
     4093  // Matrix may have been created so get rid of it
     4094  matrix_->releasePackedMatrix();
    40934095  // get matrix data pointers
    40944096  const int * column = copy.getIndices();
     
    45904592  int saveQuadraticActivated = objective_->activated();
    45914593  objective_->setActivated(0);
     4594  ClpObjective * saveObjective = objective_;
    45924595  CoinAssert (ifValuesPass>=0&&ifValuesPass<3);
    45934596  /*  Note use of "down casting".  The only class the user sees is ClpSimplex.
     
    46024605  */
    46034606  int returnCode = ((ClpSimplexDual *) this)->dual(ifValuesPass, startFinishOptions);
     4607  //int lastAlgorithm = -1;
    46044608  if ((specialOptions_&2048)!=0&&problemStatus_==10&&!numberPrimalInfeasibilities_
    46054609      &&sumDualInfeasibilities_<1000.0*dualTolerance_&&perturbation_>=100)
     
    46074611  if (problemStatus_==10) {
    46084612    //printf("Cleaning up with primal\n");
     4613    //lastAlgorithm=1;
    46094614    int savePerturbation = perturbation_;
    46104615    int saveLog = handler_->logLevel();
     
    46204625    // check which algorithms allowed
    46214626    int dummy;
    4622     if (problemStatus_==10)
     4627    if (problemStatus_==10&&saveObjective==objective_)
    46234628      startFinishOptions |= 2;
    46244629    if ((matrix_->generalExpanded(this,4,dummy)&1)!=0)
     
    46264631    else
    46274632      returnCode = ((ClpSimplexDual *) this)->dual(0,startFinishOptions);
     4633    if (saveObjective != objective_) {
     4634      // We changed objective to see if infeasible
     4635      delete objective_;
     4636      objective_=saveObjective;
     4637      if (!problemStatus_) {
     4638        // carry on
     4639        returnCode = ((ClpSimplexPrimal *) this)->primal(1,startFinishOptions);
     4640      }
     4641    }
    46284642    if (problemStatus_==3&&numberIterations_<saveMax) {
    46294643      if (handler_->logLevel()==63)
     
    46864700  //factorization_->pivotTolerance(savedPivotTolerance);
    46874701  onStopped(); // set secondary status if stopped
     4702  //if (problemStatus_==1&&lastAlgorithm==1)
     4703  //returnCode=10; // so will do primal after postsolve
    46884704  return returnCode;
    46894705}
     
    47094725  */
    47104726  int returnCode = ((ClpSimplexPrimal *) this)->primal(ifValuesPass,startFinishOptions);
     4727  //int lastAlgorithm=1;
    47114728  if (problemStatus_==10) {
     4729    //lastAlgorithm=-1;
    47124730    //printf("Cleaning up with dual\n");
    47134731    int savePerturbation = perturbation_;
     
    47344752  //factorization_->pivotTolerance(savedPivotTolerance);
    47354753  onStopped(); // set secondary status if stopped
     4754  //if (problemStatus_==1&&lastAlgorithm==1)
     4755  //returnCode=10; // so will do primal after postsolve
    47364756  return returnCode;
    47374757}
  • trunk/Clp/src/ClpSimplexDual.cpp

    r1048 r1111  
    101101#include "ClpDualRowDantzig.hpp"
    102102#include "ClpMessage.hpp"
     103#include "ClpLinearObjective.hpp"
    103104#include <cfloat>
    104105#include <cassert>
     
    15221523              // unless primal feasible!!!!
    15231524              //printf("%d %g %d %g\n",numberPrimalInfeasibilities_,sumPrimalInfeasibilities_,
    1524               //     numberDualInfeasibilities_,sumDualInfeasibilities_);
    1525               if (numberDualInfeasibilities_)
     1525              //   numberDualInfeasibilities_,sumDualInfeasibilities_);
     1526              if (sumPrimalInfeasibilities_<1.0e-3||sumDualInfeasibilities_>1.0e-6||(specialOptions_&1024)!=0) {
    15261527                problemStatus_=10;
     1528                // Get rid of objective
     1529                objective_ = new ClpLinearObjective(NULL,numberColumns_);
     1530              }
    15271531              rowArray_[0]->clear();
    15281532              columnArray_[0]->clear();
  • trunk/Clp/src/ClpSimplexNonlinear.cpp

    r1096 r1111  
    33643364  // Both these should be modified depending on problem
    33653365  // Possibly start with big bounds
    3366   double penalties[]={1.0e-3,1.0e7,1.0e9};
    3367   //double bounds[] = {1.0e-2,1.0,COIN_DBL_MAX};
     3366  //double penalties[]={1.0e-3,1.0e7,1.0e9};
     3367  double penalties[]={1.0e7,1.0e8,1.0e9};
    33683368  double bounds[] = {1.0e-2,1.0e2,COIN_DBL_MAX};
    33693369  // see how many extra we need
     
    36203620  bool zeroTargetDrop=false;
    36213621  double * gradient = new double [numberColumns_];
     3622  bool goneFeasible=false;
    36223623  // keep sum of artificials
    36233624#define KEEP_SUM 5
     
    37443745          trust[jNon] = 0.1;
    37453746    }
     3747    if (sumInfeas<0.001&&!goneFeasible) {
     3748      goneFeasible=true;
     3749      penalties[0]=1.0e-3;
     3750      penalties[1]=1.0e6;
     3751      printf("Got feasible\n");
     3752    }
    37463753    double infValue=0.0;
    37473754    double objValue=0.0;
     
    37793786      double functionValue=constraint->functionValue(this,solution);
    37803787      double dualValue = newModel.dualRowSolution()[iRow];
    3781       if (numberConstraints<50)
     3788      if (numberConstraints<-50)
    37823789        printf("For row %d current value is %g (row activity %g) , dual is %g\n",iRow,functionValue,
    37833790               newModel.primalRowSolution()[iRow],
  • trunk/Clp/src/ClpSimplexOther.cpp

    r1034 r1111  
    734734  return status;
    735735}
    736 // Creates dual of a problem
     736/* Creates dual of a problem if looks plausible
     737   (defaults will always create model)
     738   fractionRowRanges is fraction of rows allowed to have ranges
     739   fractionColumnRanges is fraction of columns allowed to have ranges
     740*/
    737741ClpSimplex *
    738 ClpSimplexOther::dualOfModel() const
     742ClpSimplexOther::dualOfModel(double fractionRowRanges,double fractionColumnRanges) const
    739743{
    740744  const ClpSimplex * model2 = (const ClpSimplex *) this;
    741745  bool changed=false;
     746  int numberChanged=0;
    742747  int iColumn;
    743748  // check if we need to change bounds to rows
     
    746751        columnLower_[iColumn]>-1.0e20) {
    747752      changed=true;
    748       break;
    749     }
     753      numberChanged++;
     754    }
     755  }
     756  int iRow;
     757  int numberExtraRows=0;
     758  if (numberChanged<=fractionColumnRanges*numberColumns_) {
     759    for (iRow=0;iRow<numberRows_;iRow++) {
     760      if (rowLower_[iRow]>-1.0e20&&
     761          rowUpper_[iRow]<1.0e20) {
     762        if (rowUpper_[iRow]!=rowLower_[iRow])
     763          numberExtraRows++;
     764      }
     765    }
     766    if (numberExtraRows>fractionRowRanges*numberRows_)
     767      return NULL;
     768  } else {
     769    return NULL;
    750770  }
    751771  if (changed) {
     
    784804  // get transpose
    785805  CoinPackedMatrix rowCopy = *matrix;
    786   int iRow;
    787   int numberExtraRows=0;
    788   for (iRow=0;iRow<numberRows;iRow++) {
    789     if (rowLower[iRow]>-1.0e20&&
    790         rowUpper[iRow]<1.0e20) {
    791       if (rowUpper[iRow]!=rowLower[iRow])
    792          numberExtraRows++;
    793     }
    794   }
    795806  const int * row = matrix->getIndices();
    796807  const int * columnLength = matrix->getVectorLengths();
     
    897908                        fromColumnsLower,fromColumnsUpper);
    898909  modelDual->setObjectiveOffset(objOffset);
     910  modelDual->setDualBound(model2->dualBound());
     911  modelDual->setInfeasibilityCost(model2->infeasibilityCost());
     912  modelDual->setDualTolerance(model2->dualTolerance());
     913  modelDual->setPrimalTolerance(model2->primalTolerance());
     914  modelDual->setPerturbation(model2->perturbation());
     915  modelDual->setSpecialOptions(model2->specialOptions());
     916  modelDual->setMoreSpecialOptions(model2->moreSpecialOptions());
    899917  delete [] fromRowsLower;
    900918  delete [] fromRowsUpper;
  • trunk/Clp/src/ClpSimplexOther.hpp

    r1034 r1111  
    163163    /// Read a basis from the given filename
    164164    int readBasis(const char *filename);
    165   /// Creates dual of a problem
    166   ClpSimplex * dualOfModel() const;
     165  /** Creates dual of a problem if looks plausible
     166      (defaults will always create model)
     167      fractionRowRanges is fraction of rows allowed to have ranges
     168      fractionColumnRanges is fraction of columns allowed to have ranges
     169  */
     170  ClpSimplex * dualOfModel(double fractionRowRanges=1.0,double fractionColumnRanges=1.0) const;
    167171  /** Restores solution from dualized problem
    168172      non-zero return code indicates minor problems
  • trunk/Clp/src/ClpSimplexPrimal.cpp

    r1095 r1111  
    13441344    }
    13451345  }
     1346  if (problemStatus_==0) {
     1347    double objVal = nonLinearCost_->feasibleCost();
     1348    double tol = 1.0e-10*CoinMax(fabs(objVal),fabs(objectiveValue_))+1.0e-8;
     1349    if (fabs(objVal-objectiveValue_)>tol) {
     1350#ifdef COIN_DEVELOP
     1351      if (handler_->logLevel()>0)
     1352        printf("nonLinearCost has feasible obj of %g, objectiveValue_ is %g\n",
     1353               objVal,objectiveValue_);
     1354#endif
     1355      objectiveValue_ = objVal;
     1356    }
     1357  }
    13461358  // save extra stuff
    13471359  matrix_->generalExpanded(this,5,dummy);
  • trunk/Clp/src/ClpSolve.cpp

    r1106 r1111  
    643643    // network - switch off stuff
    644644    doIdiot=0;
    645     doSprint=0;
     645    if (doSprint<0)
     646      doSprint=0;
    646647  }
    647648#else
     
    649650    // network - switch off stuff
    650651    doIdiot=0;
    651     doSprint=0;
     652    //doSprint=0;
    652653  }
    653654#endif
     
    848849        if(numberRows*10>numberColumns||numberColumns<6000
    849850           ||(numberRows*20>numberColumns&&!plusMinus))
    850           maxSprintPass=5;
     851          maxSprintPass=10;
    851852      }
    852853    } else if (doSprint==0) {
     
    13301331    int originalNumberColumns = model2->numberColumns();
    13311332    int numberRows = model2->numberRows();
     1333    ClpSimplex * originalModel2 = model2;
    13321334   
    13331335    // We will need arrays to choose variables.  These are too big but ..
     
    13841386    model2->clpMatrix()->times(1.0,columnSolution,rowSolution);
    13851387    // See if we can adjust using costed slacks
    1386     double penalty=CoinMin(infeasibilityCost_,1.0e10)*optimizationDirection_;
     1388    double penalty=CoinMax(1.0e5,CoinMin(infeasibilityCost_*0.01,1.0e10))*optimizationDirection_;
    13871389    const double * lower = model2->rowLower();
    13881390    const double * upper = model2->rowUpper();
     
    14261428    delete [] negSlack;
    14271429    delete [] posSlack;
    1428     int * addStarts = new int [numberRows+1];
    1429     int * addRow = new int[numberRows];
    1430     double * addElement = new double[numberRows];
     1430    int nRow=numberRows;
     1431    bool network=false;
     1432    if (dynamic_cast< ClpNetworkMatrix*>(matrix_)) {
     1433      network=true;
     1434      nRow *= 2;
     1435    }
     1436    int * addStarts = new int [nRow+1];
     1437    int * addRow = new int[nRow];
     1438    double * addElement = new double[nRow];
    14311439    addStarts[0]=0;
    14321440    int numberArtificials=0;
     1441    int numberAdd=0;
    14331442    double * addCost = new double [numberRows];
    14341443    for (iRow=0;iRow<numberRows;iRow++) {
    14351444      if (lower[iRow]>rowSolution[iRow]+1.0e-8) {
    1436         addRow[numberArtificials]=iRow;
    1437         addElement[numberArtificials]=1.0;
     1445        addRow[numberAdd]=iRow;
     1446        addElement[numberAdd++]=1.0;
     1447        if (network) {
     1448          addRow[numberAdd]=numberRows;
     1449          addElement[numberAdd++]=-1.0;
     1450        }
    14381451        addCost[numberArtificials]=penalty;
    14391452        numberArtificials++;
    1440         addStarts[numberArtificials]=numberArtificials;
     1453        addStarts[numberArtificials]=numberAdd;
    14411454      } else if (upper[iRow]<rowSolution[iRow]-1.0e-8) {
    1442         addRow[numberArtificials]=iRow;
    1443         addElement[numberArtificials]=-1.0;
     1455        addRow[numberAdd]=iRow;
     1456        addElement[numberAdd++]=-1.0;
     1457        if (network) {
     1458          addRow[numberAdd]=numberRows;
     1459          addElement[numberAdd++]=1.0;
     1460        }
    14441461        addCost[numberArtificials]=penalty;
    14451462        numberArtificials++;
    1446         addStarts[numberArtificials]=numberArtificials;
     1463        addStarts[numberArtificials]=numberAdd;
    14471464      }
    14481465    }
     
    14501467      // need copy so as not to disturb original
    14511468      model2 = new ClpSimplex(*model2);
    1452     }
    1453     model2->addColumns(numberArtificials,NULL,NULL,addCost,
    1454                        addStarts,addRow,addElement);
     1469      if (network) {
     1470        // network - add a null row
     1471        model2->addRow(0,NULL,NULL,-COIN_DBL_MAX,COIN_DBL_MAX);
     1472        numberRows++;
     1473      }
     1474      model2->addColumns(numberArtificials,NULL,NULL,addCost,
     1475                         addStarts,addRow,addElement);
     1476    }
    14551477    delete [] addStarts;
    14561478    delete [] addRow;
     
    15671589    model2->getDblParam(ClpObjOffset,originalOffset);
    15681590    int totalIterations=0;
     1591    double lastSumArtificials=COIN_DBL_MAX;
     1592    int originalMaxSprintPass=maxSprintPass;
     1593    maxSprintPass=20; // so we do that many if infeasible
    15691594    for (iPass=0;iPass<maxSprintPass;iPass++) {
    15701595      //printf("Bug until submodel new version\n");
     
    16651690      CoinMemcpyN(small.primalRowSolution(),
    16661691             numberRows,model2->primalRowSolution());
     1692      double sumArtificials = 0.0;
     1693      for (i=0;i<numberArtificials;i++)
     1694        sumArtificials += fullSolution[i + originalNumberColumns];
     1695      if (sumArtificials&&iPass>5&&sumArtificials>=lastSumArtificials) {
     1696        // increase costs
     1697        double * cost = model2->objective()+originalNumberColumns;
     1698        double newCost = CoinMin(1.0e10,cost[0]*1.5);
     1699        for (i=0;i<numberArtificials;i++)
     1700          cost[i]=newCost;
     1701      }
     1702      lastSumArtificials = sumArtificials;
    16671703      // get reduced cost for large problem
    16681704      double * djs = model2->dualColumnSolution();
     
    16911727        <<numberNegative
    16921728        <<CoinMessageEol;
    1693       if ((small.objectiveValue()*optimizationDirection_>lastObjective-1.0e-7&&iPass>5)||
     1729      if (sumArtificials<1.0e-8&&originalMaxSprintPass>=0) {
     1730        maxSprintPass = iPass+originalMaxSprintPass;
     1731        originalMaxSprintPass=-1;
     1732      }
     1733      if (iPass>20)
     1734        sumArtificials=0.0;
     1735      if ((small.objectiveValue()*optimizationDirection_>lastObjective-1.0e-7&&iPass>5&&sumArtificials<1.0e-8)||
    16941736          (!small.numberIterations()&&iPass)||
    16951737          iPass==maxSprintPass-1||small.status()==3) {
     
    17331775      sort[i] = i + originalNumberColumns;
    17341776    model2->deleteColumns(numberArtificials,sort);
     1777    if (network) {
     1778      int iRow=numberRows-1;
     1779      model2->deleteRows(1,&iRow);
     1780    }
    17351781    delete [] weight;
    17361782    delete [] sort;
     
    17531799    }
    17541800    model2->primal(1);
     1801    if (model2!=originalModel2)
     1802      originalModel2->moveInfo(*model2);
    17551803    model2->setPerturbation(savePerturbation);
    17561804    time2 = CoinCpuTime();
Note: See TracChangeset for help on using the changeset viewer.