Changeset 2074 for trunk


Ignore:
Timestamp:
Dec 10, 2014 4:43:54 AM (5 years ago)
Author:
forrest
Message:

changes to try and make more stable

Location:
trunk/Clp/src
Files:
24 edited

Legend:

Unmodified
Added
Removed
  • trunk/Clp/src/AbcCommon.hpp

    r1910 r2074  
    3838#ifndef ABC_INHERIT
    3939#define ABC_INHERIT
     40#ifndef CLP_INHERIT_MODE
     41#define CLP_INHERIT_MODE 2
     42#endif
    4043#endif
    4144#elif CLP_HAS_ABC==3
  • trunk/Clp/src/AbcMatrix.cpp

    r2024 r2074  
    31893189   int maximumRows=model_->maximumAbcNumberRows();
    31903190   int number=updateForTableauRow.getNumElements();
    3191    // was assert (number);
     3191#ifdef GCC_4_9
     3192   assert (number);
     3193#else
    31923194   if (!number) {
    31933195     printf("Null tableau row!\n");
    31943196   }
     3197#endif
    31953198   bool useRowCopy = (gotRowCopy()&&(number<<2)<maximumRows);
    31963199   //useRowCopy=false;
  • trunk/Clp/src/AbcSimplex.cpp

    r2032 r2074  
    14961496      maximumPivots < 2) {
    14971497    // If dense then increase
    1498     if (maximumPivots > 100 && numberDense > 1.5 * maximumPivots) {
     1498    if (maximumPivots > 100 && numberDense > 1.5 * maximumPivots && false) {
    14991499      abcFactorization_->maximumPivots(numberDense);
    15001500    }
  • trunk/Clp/src/AbcSimplexPrimal.cpp

    r2030 r2074  
    37943794        // check update
    37953795        CoinIndexedVector * tempVector[4];
    3796         memcpy(tempVector,vector+4*iBest,4*sizeof(CoinIndexedVector *));
     3796        for (int i=0;i<4;i++)
     3797          tempVector[i] = vector[4*iBest+i];
    37973798        returnCode = cilk_spawn doFTUpdate(tempVector);
    37983799        bestUpdate=tempVector[0];
  • trunk/Clp/src/CbcOrClpParam.cpp

    r2070 r2074  
    17861786     parameters[numberParameters-1].append("so!low_halim");
    17871787     parameters[numberParameters-1].append("lots");
    1788 #ifdef ABC_INHERIT
     1788#ifdef CLP_INHERIT_MODE
    17891789     parameters[numberParameters-1].append("dual");
    17901790     parameters[numberParameters-1].append("dw");
  • trunk/Clp/src/ClpMain.cpp

    r2030 r2074  
    1717void ClpMain0(AbcSimplex * models);
    1818int ClpMain1(int argc, const char *argv[],AbcSimplex * model);
    19 #endif
     19#endif 
    2020//#define CILK_TEST
    2121#ifdef CILK_TEST
    2222static void cilkTest();
     23#endif
     24//#define LAPACK_TEST
     25//#define CLP_USE_OPENBLAS 1
     26#if CLP_USE_OPENBLAS
     27extern "C"
     28{
     29  void openblas_set_num_threads(int num_threads);
     30}
     31#endif
     32#ifdef LAPACK_TEST
     33#include "/include/lapacke.h"
     34#ifndef COIN_FACTORIZATION_DENSE_CODE
     35#define COIN_FACTORIZATION_DENSE_CODE 1
     36#endif
     37#if COIN_FACTORIZATION_DENSE_CODE==1
     38// using simple lapack interface
     39extern "C"
     40{
     41  void openblas_set_num_threads(int num_threads);
     42#if 0
     43  /** LAPACK Fortran subroutine DGETRF. */
     44  void LAPACK_dgetrf(int * m, int *n,
     45                               double *A, int *ldA,
     46                               int * ipiv, int *info);
     47  /** LAPACK Fortran subroutine DGETRS. */
     48  void LAPACK_dgetrs(char *trans, int *n,
     49                               int *nrhs, const double *A, int *ldA,
     50                     int * ipiv, double *B, int *ldB, int *info);
     51  //    LAPACK_dgetrf(&N, &N, m, &LDA,ipiv, &info);
     52#endif
     53}
     54int test_lapack(int n)
     55{
     56  int* ipiv;
     57  int info;
     58  int i, j;
     59  double * m, *x, *y;
     60 
     61  int LDB,LDA, N, NRHS;
     62  char transp = 'N';
     63 
     64 
     65  m=(double*)malloc(sizeof(double)*n*n);
     66  x=(double*)malloc(sizeof(double)*n);
     67  y=(double*)malloc(sizeof(double)*n);
     68  ipiv=(int*)malloc(sizeof(int)*n);
     69 
     70  for (i=0; i<n; ++i) {
     71    x[i]=1.0;
     72    for (j=0; j<n; ++j) {
     73      m[i*n+j]=(rand()%100+1)/10.0;
     74      //      printf("m[%d,%d]=%lf\n",i,j, m[i*n+j]);
     75    }
     76  }
     77 
     78  /* test cblas.h */
     79  //cblas_dgemv(CblasColMajor, CblasNoTrans, n, n, 1.0, m, n,
     80  //          x, 1, 0.0, y, 1);
     81 
     82  //  for (i=0; i<n; ++i)  printf("x[%d]=%lf\n",i, x[i]);
     83  //for (i=0; i<n; ++i)  printf("y[%d]=%lf\n",i, y[i]);
     84 
     85        LDB=n;
     86        LDA=n;
     87        N=n;
     88        NRHS=1;
     89        info=0;
     90       
     91        LAPACK_dgetrf(&N, &N, m, &LDA,ipiv, &info);
     92       
     93if (info != 0) fprintf(stderr, "dgetrf failure with error %d\n", info);
     94 
     95  LAPACK_dgetrs(&transp, &N, &NRHS, m, &LDA, ipiv, y, &LDB, &info);
     96 
     97  if (info != 0) fprintf(stderr, "failure with error %d\n", info);
     98  //  for (i=0; i<n; ++i) printf("%lf\n", y[i]);
     99
     100  free(m);
     101  free(x);
     102  free(y);
     103  free(ipiv);
     104  return 0;
     105}
     106#elif COIN_FACTORIZATION_DENSE_CODE==2
     107// C interface
     108enum CBLAS_ORDER {CblasRowMajor=101, CblasColMajor=102};
     109enum CBLAS_TRANSPOSE {CblasNoTrans=111, CblasTrans=112};
     110extern "C"
     111{
     112int clapack_dgetrf ( const enum CBLAS_ORDER Order, const int M, const int N, double *A, const int lda, int *ipiv );
     113int clapack_dgetrs ( const enum CBLAS_ORDER Order,
     114                       const enum CBLAS_TRANSPOSE Trans,
     115                       const int N, const int NRHS,
     116                       const double *A, const int lda, const int *ipiv, double *B,
     117                       const int ldb );
     118}
     119int test_lapack(int n)
     120{
     121  int* ipiv;
     122  int info;
     123  int i, j;
     124  double * m, *x, *y;
     125 
     126  int LDB,LDA, N, NRHS;
     127  char transp = 'N';
     128 
     129 
     130  m=(double*)malloc(sizeof(double)*n*n);
     131  x=(double*)malloc(sizeof(double)*n);
     132  y=(double*)malloc(sizeof(double)*n);
     133  ipiv=(int*)malloc(sizeof(int)*n);
     134 
     135  for (i=0; i<n; ++i) {
     136    x[i]=1.0;
     137    for (j=0; j<n; ++j) {
     138      m[i*n+j]=(rand()%100+1)/10.0;
     139      //      printf("m[%d,%d]=%lf\n",i,j, m[i*n+j]);
     140    }
     141  }
     142 
     143  /* test cblas.h */
     144  //cblas_dgemv(CblasColMajor, CblasNoTrans, n, n, 1.0, m, n,
     145  //          x, 1, 0.0, y, 1);
     146 
     147  //  for (i=0; i<n; ++i)  printf("x[%d]=%lf\n",i, x[i]);
     148  //for (i=0; i<n; ++i)  printf("y[%d]=%lf\n",i, y[i]);
     149 
     150        LDB=n;
     151        LDA=n;
     152        N=n;
     153        NRHS=1;
     154        info=clapack_dgetrf ( CblasColMajor,n,n,m,n,ipiv );
     155       
     156if (info != 0) fprintf(stderr, "dgetrf failure with error %d\n", info);
     157 
     158  clapack_dgetrs ( CblasColMajor,CblasNoTrans,n,1,m,n,ipiv,y,n);
     159 
     160  if (info != 0) fprintf(stderr, "failure with error %d\n", info);
     161  //  for (i=0; i<n; ++i) printf("%lf\n", y[i]);
     162
     163  free(m);
     164  free(x);
     165  free(y);
     166  free(ipiv);
     167  return 0;
     168}
     169#endif
    23170#endif
    24171int
     
    28175main (int argc, const char *argv[])
    29176{
    30 #ifdef CILK_TEST
     177#ifdef CILK_TEST 
    31178  cilkTest();
     179#endif 
     180#if CLP_USE_OPENBLAS
     181  openblas_set_num_threads(CLP_USE_OPENBLAS);
     182#endif 
     183#ifdef LAPACK_TEST
     184  //void openblas_set_num_threads(int num_threads);
     185  openblas_set_num_threads(1);
     186  if(argc<2){
     187    printf("Error - need size of matrix for lapack test\n");
     188    return 1;
     189  }
     190  int n=atoi(argv[1]);
     191  printf("n=%d\n",n);
     192  if(argc>2){
     193    int nThreads=atoi(argv[2]);
     194    printf("Using %d threads\n",nThreads);
     195    openblas_set_num_threads(nThreads);
     196  }
     197  test_lapack(n);
     198  return 0;
    32199#endif
    33200#ifndef ABC_INHERIT
  • trunk/Clp/src/ClpModel.hpp

    r2046 r2074  
    10411041         2097152 - zero costs!
    10421042         4194304 - don't scale integer variables
     1043         8388608 - Idiot when not really sure about it
    10431044         NOTE - many applications can call Clp but there may be some short cuts
    10441045                which are taken which are not guaranteed safe from all applications.
  • trunk/Clp/src/ClpPackedMatrix.cpp

    r1972 r2074  
    12131213               }
    12141214#endif
    1215                //#define COIN_SPARSE_MATRIX 1
    1216 #if COIN_SPARSE_MATRIX
    1217                assert (!y->getNumElements());
     1215#define COIN_SPARSE_MATRIX 2
     1216               int numberCovered=0;
     1217               int numberColumns = matrix_->getNumCols();
     1218               bool sparse=true;
     1219               for (int i = 0; i < numberInRowArray; i++) {
     1220                 int iRow = whichRow[i];
     1221                 numberCovered += rowStart[iRow+1] - rowStart[iRow];
     1222                 // ? exact ratio
     1223                 if (numberCovered>numberColumns) {
     1224                   sparse=false;
     1225                   break;
     1226                 }
     1227               }
     1228               if (sparse) {
     1229                 assert (!y->getNumElements());
    12181230#if COIN_SPARSE_MATRIX != 2
    1219               // and set up mark as char array
    1220               char * marked = reinterpret_cast<char *> (index+columnArray->capacity());
    1221               int * lookup = y->getIndices();
     1231                // and set up mark as char array
     1232                char * marked = reinterpret_cast<char *> (index+columnArray->capacity());
     1233                int * lookup = y->getIndices();
    12221234#ifndef NDEBUG
    1223               //int numberColumns = matrix_->getNumCols();
    1224               //for (int i=0;i<numberColumns;i++)
    1225               //assert(!marked[i]);
    1226 #endif
    1227               numberNonZero=gutsOfTransposeTimesByRowGE3a(rowArray,index,array,
    1228                                            lookup,marked,zeroTolerance,scalar);
     1235                //int numberColumns = matrix_->getNumCols();
     1236                //for (int i=0;i<numberColumns;i++)
     1237                //assert(!marked[i]);
     1238#endif
     1239                numberNonZero=gutsOfTransposeTimesByRowGE3a(rowArray,index,array,
     1240                                                             lookup,marked,zeroTolerance,scalar);
    12291241#else
    1230                double  * array2 = y->denseVector();
    1231                numberNonZero=gutsOfTransposeTimesByRowGE3(rowArray,index,array,
    1232                                            array2,zeroTolerance,scalar);
    1233 #endif
    1234 #else
    1235                int numberColumns = matrix_->getNumCols();
    1236                numberNonZero = gutsOfTransposeTimesByRowGEK(rowArray, index, array,
    1237                                numberColumns, zeroTolerance, scalar);
    1238 #endif
     1242                 double  * array2 = y->denseVector();
     1243                 numberNonZero=gutsOfTransposeTimesByRowGE3(rowArray,index,array,
     1244                                                            array2,zeroTolerance,scalar);
     1245#endif
     1246               } else {
     1247                 numberNonZero = gutsOfTransposeTimesByRowGEK(rowArray, index, array,
     1248                                                              numberColumns, zeroTolerance, scalar);
     1249               }
    12391250               columnArray->setNumElements(numberNonZero);
    12401251          } else {
  • trunk/Clp/src/ClpPresolve.cpp

    r2070 r2074  
    5959     nrows_(0),
    6060     nelems_(0),
    61 #ifdef ABC_INHERIT
     61#ifdef CLP_INHERIT_MODE
    6262     numberPasses_(20),
    6363#else
     
    893893     CoinMessages messages = CoinMessage(prob->messages().language());
    894894     paction_ = 0;
    895      prob->maxSubstLevel_ = 3 ;
     895     prob->maxSubstLevel_ = CoinMax(3,prob->maxSubstLevel_) ;
    896896#ifndef PRESOLVE_DETAIL
    897897     if (prob->tuning_) {
     
    956956          if (doIntersection())
    957957            prob->presolveOptions_ |= 0x10;
    958 
     958          // zero small elements in aggregation
     959          prob->presolveOptions_ |= zeroSmall()*0x20000;
    959960          // some things are expensive so just do once (normally)
    960961
     
    995996               printProgress('C',0);
    996997          }
    997 #ifdef ABC_INHERIT
    998998          if (doTwoxTwo()) {
    999999            possibleSkip;
    10001000            paction_ = twoxtwo_action::presolve(prob, paction_);
    10011001          }
    1002 #endif
    10031002          if (duprow) {
    10041003            possibleSkip;
     
    10261025          int lastDropped = 0;
    10271026          prob->pass_ = 0;
    1028 #ifdef ABC_INHERIT
     1027#if CLP_INHERIT_MODE>1
    10291028          int numberRowsStart=nrows_-prob->countEmptyRows();
    10301029          int numberColumnsStart=ncols_-prob->countEmptyCols();
     
    13381337               if (paction_ == paction0 || stopLoop)
    13391338                    break;
    1340 #ifdef ABC_INHERIT
     1339#if CLP_INHERIT_MODE>1
    13411340               // see whether to stop anyway
    13421341               int numberRowsNow=nrows_-prob->countEmptyRows();
     
    14371436                         double coeff = colels[k];
    14381437                         k = link[k];
     1438                         assert (k!=NO_LINK||i==nx-1);
    14391439                         rsol[row] += solutionValue * coeff;
    14401440                    }
     
    16351635
    16361636{
    1637      bulk0_ = static_cast<CoinBigIndex> (bulkRatio_ * nelems_in);
     1637     bulk0_ = static_cast<CoinBigIndex> (bulkRatio_ *
     1638                                         CoinMax(nelems_in,nelems_));
    16381639     hrow_  = new int   [bulk0_];
    16391640     colels_ = new double[bulk0_];
     
    21972198          double ratio = 2.0;
    21982199          if (substitution_ > 3)
    2199                ratio = substitution_;
     2200               ratio = sqrt((substitution_-3)+5.0);
    22002201          else if (substitution_ == 2)
    22012202               ratio = 1.5;
  • trunk/Clp/src/ClpPresolve.hpp

    r1972 r2074  
    184184          else presolveActions_ |= 4096;
    185185     }
     186     /** How much we want to zero small values from aggregation - ratio
     187         0 - 1.0e-12, 1 1.0e-11, 2 1.0e-10, 3 1.0e-9 */
     188     inline int zeroSmall() const {
     189          return (presolveActions_&(8192|16384))>>13;
     190     }
     191     inline void setZeroSmall(int value) {
     192         presolveActions_  &= ~(8192|16384);
     193         presolveActions_ |= value<<13;
     194     }
    186195     /// Set whole group
    187196     inline int presolveActions() const {
  • trunk/Clp/src/ClpPrimalColumnSteepest.cpp

    r2068 r2074  
    28722872{
    28732873     model_ = model;
     2874     if (mode==6) {
     2875       // If incoming weight is 1.0 then return else as 5
     2876       assert (weights_);
     2877       int sequenceIn = model_->sequenceIn();
     2878       assert (sequenceIn>=0&&sequenceIn<model_->numberRows()+model_->numberColumns());
     2879       if (weights_[sequenceIn]==(mode_!=1) ? 1.0 : 1.0+ADD_ONE)
     2880         return;
     2881       else
     2882         mode=5;
     2883     }
    28742884     if (mode_ == 4 || mode_ == 5) {
    28752885          if (mode == 1 && !weights_)
  • trunk/Clp/src/ClpSimplex.cpp

    r2072 r2074  
    544544     double * save = NULL;
    545545     double oldValue = 0.0;
     546     double oldLargestPrimalError=largestPrimalError_;
     547     double oldLargestDualError=largestDualError_;
    546548     if (valuesPass) {
    547549          assert(algorithm_ > 0); // only primal at present
     
    756758               factorization_->zeroTolerance(1.0e-18);
    757759     }
     760     int returnCode=0;
     761     if (numberIterations_ && (forceFactorization_ > 2 || forceFactorization_<0 || factorization_->pivotTolerance()<0.9899999999)) {
     762       if ((largestPrimalError_>1.0e3&&
     763           oldLargestPrimalError*1.0e2<largestPrimalError_)||
     764           (largestDualError_>1.0e3&&
     765            oldLargestDualError*1.0e2<largestDualError_)) {
     766         double pivotTolerance = factorization_->pivotTolerance();
     767         double factor=(largestPrimalError_>1.0e10||largestDualError_>1.0e10)
     768           ? 2.0 : 1.2;
     769         if (pivotTolerance<0.1)
     770           factorization_->pivotTolerance(0.1);
     771         else if (pivotTolerance<0.98999999)
     772           factorization_->pivotTolerance(CoinMin(0.99,pivotTolerance*factor));
     773#ifdef CLP_USEFUL_PRINTOUT
     774         if (pivotTolerance<0.9899999) {
     775           printf("Changing pivot tolerance from %g to %g and backtracking\n",
     776                  pivotTolerance,factorization_->pivotTolerance());
     777         }
     778         printf("because old,new primal error %g,%g - dual %g,%g pivot_tol %g\n",
     779                oldLargestPrimalError,largestPrimalError_,
     780                oldLargestDualError,largestDualError_,
     781                pivotTolerance);
     782#endif
     783         if (pivotTolerance<0.9899999) {
     784           largestPrimalError_=0.0;
     785           largestDualError_=0.0;
     786           returnCode=1;
     787         }
     788       }
     789     }
    758790     // Switch off false values pass indicator
    759791     if (!valuesPass && algorithm_ > 0)
     
    764796            numberPrimalInfeasibilities_,sumPrimalInfeasibilities_,sumOfRelaxedPrimalInfeasibilities_,
    765797            numberDualInfeasibilities_,sumDualInfeasibilities_,sumOfRelaxedDualInfeasibilities_);
    766      return 0;
     798     return returnCode;
    767799}
    768800void
     
    852884     CoinIndexedVector * thisVector = arrayVector;
    853885     CoinIndexedVector * lastVector = previousVector;
     886#if 0
     887     static double * xsave=NULL;
     888     {
     889       double * xx = thisVector->denseVector();
     890       double largest=0.0;
     891       int iLargest=-1;
     892       for (int i=0;i<numberRows_;i++) {
     893         if (fabs(xx[i])>largest) {
     894           largest=fabs(xx[i]);
     895           iLargest=i;
     896         }
     897       }
     898       printf("largest incoming rhs %g on row %d\n",largest,iLargest);
     899     }
     900     if (numberIterations_<-40722) {
     901       double * xx = thisVector->denseVector();
     902       if (xsave) {
     903         double * sol = xsave+numberRows_;
     904         double largest=0.0;
     905         int iLargest=-1;
     906         for (int i=0;i<numberRows_;i++) {
     907           if (fabs(xx[i]-xsave[i])>largest) {
     908             largest=fabs(xx[i]-xsave[i]);
     909             iLargest=i;
     910           }
     911         }
     912         printf("error %g on row %d\n",largest,iLargest);
     913         largest=0.0;
     914         iLargest=-1;
     915         for (int i=0;i<numberColumns_;i++) {
     916           if (fabs(solution_[i]-sol[i])>largest) {
     917             largest=fabs(solution_[i]-sol[i]);
     918             iLargest=i;
     919           }
     920         }
     921         printf("error %g on col %d\n",largest,iLargest);
     922       } else {
     923         xsave=new double[2*numberRows_+numberColumns_];
     924       }
     925       memcpy(xsave,xx,numberRows_*sizeof(double));
     926       memcpy(xsave+numberRows_,solution_,(numberRows_+numberColumns_)*sizeof(double));
     927     }
     928#endif
    854929     if (number)
    855930          factorization_->updateColumn(workSpace, thisVector);
     
    18311906                    << status
    18321907                    << CoinMessageEol;
     1908#ifdef CLP_USEFUL_PRINTOUT
     1909          printf("Basis singular - pivot tolerance %g\n",
     1910                 factorization_->pivotTolerance());
     1911#endif
    18331912          return -1;
    18341913     } else if (!solveType) {
     
    20842163               maximumPivots < 2) {
    20852164          // If dense then increase
    2086           if (maximumPivots > 100 && numberDense > 1.5 * maximumPivots) {
     2165          if (maximumPivots > 100 && numberDense > 1.5 * maximumPivots
     2166              && false) {
    20872167               factorization_->maximumPivots(numberDense);
    20882168               dualRowPivot_->maximumPivotsChanged();
     
    21082188     } else if (numberIterations_ > 1000 + 10 * (numberRows_ + (numberColumns_ >> 2)) && matrix_->type()<15) {
    21092189          double random = randomNumberGenerator_.randomDouble();
     2190          while (random<0.45)
     2191            random *= 2.0;
    21102192          int maxNumber = (forceFactorization_ < 0) ? maximumPivots : CoinMin(forceFactorization_, maximumPivots);
    21112193          if (factorization_->pivots() >= random * maxNumber) {
     
    86058687               // Switch off dense (unless special option set)
    86068688               if ((specialOptions_ & 8) == 0)
    8607                     factorization_->setDenseThreshold(0);
     8689                  factorization_->setDenseThreshold(-saveThreshold);
    86088690          }
    86098691          // If values pass then perturb (otherwise may be optimal so leave a bit)
  • trunk/Clp/src/ClpSimplex.hpp

    r2071 r2074  
    3434 */
    3535#ifdef ABC_INHERIT
     36#ifndef CLP_INHERIT_MODE
     37#define CLP_INHERIT_MODE 1
     38#endif
    3639#ifndef ABC_CLP_DEFAULTS
    3740#define ABC_CLP_DEFAULTS 0
  • trunk/Clp/src/ClpSimplexDual.cpp

    r2070 r2074  
    15211521                              returnCode = -2;
    15221522                              moreSpecialOptions_ |= 16;
     1523                              double pivotTolerance = factorization_->pivotTolerance();
     1524                              if (pivotTolerance<0.4&&factorization_->pivots()<100) {
     1525                                factorization_->pivotTolerance(1.05*pivotTolerance);
     1526#ifdef CLP_USEFUL_PRINTOUT
     1527                                printf("Changing pivot tolerance from %g to %g as ftran/btran error %g/%g\n",
     1528                                       pivotTolerance,factorization_->pivotTolerance(),
     1529                                       alpha_,btranAlpha);
     1530#endif
     1531                              }
    15231532                              break;
    15241533                         } else {
     
    44764485     } else if (goodAccuracy()) {
    44774486          // Can reduce tolerance
    4478           double newTolerance = CoinMax(0.99 * factorization_->pivotTolerance(), saveData.pivotTolerance_);
     4487          double newTolerance = CoinMax(0.995 * factorization_->pivotTolerance(), saveData.pivotTolerance_);
    44794488          factorization_->pivotTolerance(newTolerance);
    44804489     }
     
    60676076               }
    60686077          }
     6078     }
     6079     if (largestZero>1.0*largest&&largest) {
     6080       //printf("largest zero perturbation of %g too big (nonzero %g)\n",
     6081       //     largestZero,largest);
     6082       largestZero = 0.0;
     6083       const double * obj = objective();
     6084       double test=CoinMax(1.0e-8,largest);
     6085       for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
     6086         if (!obj[iColumn]) {
     6087           double cost = cost_[iColumn];
     6088           while (fabs(cost)>test)
     6089             cost *= 0.5;
     6090           cost_[iColumn] = cost;
     6091           largestZero=CoinMax(largestZero,fabs(cost));
     6092         }
     6093       }
    60696094     }
    60706095     handler_->message(CLP_SIMPLEX_PERTURB, messages_)
  • trunk/Clp/src/ClpSimplexOther.cpp

    r2070 r2074  
    11121112                         columnActivity_[iColumn] = columnUpper_[iColumn];
    11131113                    } else {
    1114                          abort();
     1114                         setColumnStatus(iColumn, superBasic);
     1115                         columnActivity_[iColumn] = otherValue;
    11151116                    }
    11161117               }
     
    11551156                              columnActivity_[iColumn] = columnUpper_[iColumn];
    11561157                         } else {
    1157                               abort();
     1158                              setColumnStatus(iColumn, superBasic);
     1159                              columnActivity_[iColumn] = otherValue;
    11581160                         }
    11591161                    }
  • trunk/Clp/src/ClpSimplexPrimal.cpp

    r2070 r2074  
    541541               // Iterate
    542542               whileIterating(ifValuesPass ? 1 : 0);
     543               if (sequenceIn_<0&&ifValuesPass==2)
     544                 problemStatus_=3; // user wants to exit
    543545          }
    544546     }
     
    958960          }
    959961#endif
     962          if (ifValuesPass && firstFree_ <0) {
     963            largestPrimalError_=1.0e7;
     964            largestDualError_=1.0e7;
     965          }
    960966          numberThrownOut = gutsOfSolution(NULL, NULL, (firstFree_ >= 0));
    961967          double sumInfeasibility =  nonLinearCost_->sumInfeasibilities();
     
    15991605               sumPrimalInfeasibilities_ = nonLinearCost_->sumInfeasibilities();
    16001606               // say infeasible
    1601                if (numberPrimalInfeasibilities_)
     1607               if (numberPrimalInfeasibilities_ && largestPrimalError_<1.0e-1)
    16021608                    problemStatus_ = 1;
    16031609          }
     
    25002506          if (perturbation_ > 50 && perturbation_ < 55) {
    25012507               // reduce
    2502                while (perturbation_ > 50) {
    2503                     perturbation_--;
     2508               while (perturbation_ < 55) {
     2509                    perturbation_++;
    25042510                    perturbation *= 0.25;
    25052511                    bias *= 0.25;
     
    25572563     if (type == 1) {
    25582564          double tolerance = 100.0 * primalTolerance_;
     2565          tolerance = 10.0 * primalTolerance_; // try smaller
    25592566          //double multiplier = perturbation*maximumFraction;
    25602567          for (iSequence = 0; iSequence < numberRows_ + numberColumns_; iSequence++) {
     
    25752582                         value *= perturbationArray_[2*iSequence];
    25762583#endif
     2584                         if (value) {
     2585                           while (value<tolerance)
     2586                             value *= 3.0;
     2587                         }
    25772588                         if (solutionValue - lowerValue <= primalTolerance_) {
    25782589                              lower_[iSequence] -= value;
     
    26112622     } else {
    26122623          double tolerance = 100.0 * primalTolerance_;
     2624          tolerance = 10.0 * primalTolerance_; // try smaller
    26132625          for (i = 0; i < numberColumns_; i++) {
    26142626               double lowerValue = lower_[i], upperValue = upper_[i];
     
    26232635                    value *= randomNumberGenerator_.randomDouble();
    26242636                    if (savePerturbation != 50) {
    2625                          if (fabs(value) <= primalTolerance_)
    2626                               value = 0.0;
    2627                          if (lowerValue > -1.0e20 && lowerValue)
    2628                               lowerValue -= value * (CoinMax(1.0e-2, 1.0e-5 * fabs(lowerValue)));
    2629                          if (upperValue < 1.0e20 && upperValue)
    2630                               upperValue += value * (CoinMax(1.0e-2, 1.0e-5 * fabs(upperValue)));
    2631                     } else if (value) {
     2637                      if (fabs(value) <= primalTolerance_)
     2638                        value = 0.0;
     2639                    }
     2640                    if (value) {
    26322641                         double valueL = value * (CoinMax(1.0e-2, 1.0e-5 * fabs(lowerValue)));
    26332642                         // get in range
     
    26362645                              while (valueL <= tolerance)
    26372646                                   valueL *= 10.0;
    2638                          } else if (valueL > 1.0) {
     2647                         } else if (valueL > 1.0e-3) {
    26392648                              valueL *= 0.1;
    2640                               while (valueL > 1.0)
     2649                              while (valueL > 1.0e-3)
    26412650                                   valueL *= 0.1;
    26422651                         }
     
    26492658                              while (valueU <= tolerance)
    26502659                                   valueU *= 10.0;
    2651                          } else if (valueU > 1.0) {
     2660                         } else if (valueU > 1.0e-3) {
    26522661                              valueU *= 0.1;
    2653                               while (valueU > 1.0)
     2662                              while (valueU > 1.0e-3)
    26542663                                   valueU *= 0.1;
    26552664                         }
     
    30113020                    } else {
    30123021                         // take on more relaxed criterion
    3013                          if (saveDj * dualIn_ < test1 ||
     3022                         if ((saveDj * dualIn_ < test1 ||
    30143023                                   fabs(saveDj - dualIn_) > 2.0e-1 * (1.0 + fabs(dualIn_)) ||
    3015                                    fabs(dualIn_) < test2) {
     3024                              fabs(dualIn_) < test2) &&
     3025                             (fabs(saveDj)>fabs(dualIn_)
     3026                                                ||saveDj*dualIn_<1.0e-4||factorization_->pivots())) {
    30163027                              // need to reject something
    30173028                              char x = isColumn(sequenceIn_) ? 'C' : 'R';
     
    30213032                              setFlagged(sequenceIn_);
    30223033#if 1 //def FEB_TRY
     3034                              // could do conditional reset of weights to get larger djs
     3035                              primalColumnPivot_->saveWeights(this,6);
    30233036                              // Make safer?
     3037                              double oldTolerance=factorization_->pivotTolerance();
    30243038                              factorization_->saferTolerances (-0.99, -1.03);
     3039                              if (factorization_->pivotTolerance()<1.029*oldTolerance
     3040                                  &&oldTolerance<0.995
     3041                                  &&!factorization_->pivots()) {
     3042#ifdef CLP_USEFUL_PRINTOUT
     3043                                printf("Changing pivot tolerance from %g to %g and factorizing\n",
     3044                                       oldTolerance,factorization_->pivotTolerance());
     3045#endif
     3046                                clearAll();
     3047                                pivotRow_ = -1; // say no weights update
     3048                                returnCode = -4;
     3049                                if(lastGoodIteration_ + 1 == numberIterations_) {
     3050                                  // not looking wonderful - try cleaning bounds
     3051                                  // put non-basics to bounds in case tolerance moved
     3052                                  nonLinearCost_->checkInfeasibilities(0.0);
     3053                                }
     3054                                sequenceOut_ = -1;
     3055                                break;
     3056                              }
    30253057#endif
    30263058                              progress_.clearBadTimes();
  • trunk/Clp/src/ClpSolve.cpp

    r2071 r2074  
    5454#include "CoinStructuredModel.hpp"
    5555double zz_slack_value=0.0;
     56#ifdef CLP_USEFUL_PRINTOUT
     57double debugDouble[10];
     58int debugInt[24];
     59#endif
    5660
    5761#include "ClpPresolve.hpp"
     
    624628#endif
    625629    }
     630    int numberCpu=this->abcState()&15;
     631    if (numberCpu==9) {
     632      numberCpu=1;
     633#if ABC_PARALLEL==2
     634#ifndef FAKE_CILK
     635      if (number_cilk_workers>1)
     636      numberCpu=CoinMin(2*number_cilk_workers,8);
     637#endif
     638#endif
     639    } else if (numberCpu==10) {
     640      // maximum
     641      numberCpu=4;
     642    } else if (numberCpu==10) {
     643      // decide
     644      if (abcModel2->getNumElements()<5000)
     645        numberCpu=1;
     646#if ABC_PARALLEL==2
     647#ifndef FAKE_CILK
     648      else if (number_cilk_workers>1)
     649        numberCpu=CoinMin(2*number_cilk_workers,8);
     650#endif
     651#endif
     652      else
     653        numberCpu=1;
     654    } else {
     655#if ABC_PARALLEL==2
     656#ifndef FAKE_CILK
     657      char temp[3];
     658      sprintf(temp,"%d",numberCpu);
     659      __cilkrts_set_param("nworkers",temp);
     660      printf("setting cilk workers to %d\n",numberCpu);
     661      number_cilk_workers=numberCpu;
     662
     663#endif
     664#endif
     665    }
    626666    char line[200];
    627667#if ABC_PARALLEL
     
    637677#endif
    638678#endif
    639     int numberCpu=this->abcState()&15;
    640     if (numberCpu==9) {
    641       numberCpu=1;
    642 #if ABC_PARALLEL==2
    643 #ifndef FAKE_CILK
    644       if (number_cilk_workers>1)
    645       numberCpu=CoinMin(2*number_cilk_workers,8);
    646 #endif
    647 #endif
    648     } else if (numberCpu==10) {
    649       // maximum
    650       numberCpu=4;
    651     } else if (numberCpu==10) {
    652       // decide
    653       if (abcModel2->getNumElements()<5000)
    654         numberCpu=1;
    655 #if ABC_PARALLEL==2
    656 #ifndef FAKE_CILK
    657       else if (number_cilk_workers>1)
    658         numberCpu=CoinMin(2*number_cilk_workers,8);
    659 #endif
    660 #endif
    661       else
    662         numberCpu=1;
    663     } else {
    664 #if ABC_PARALLEL==2
    665 #if 0 //ndef FAKE_CILK
    666       char temp[3];
    667       sprintf(temp,"%d",numberCpu);
    668       __cilkrts_set_param("nworkers",temp);
    669 #endif
    670 #endif
    671     }
    672679    abcModel2->setParallelMode(numberCpu-1);
    673680#endif
     
    843850     ClpMatrixBase * saveMatrix = NULL;
    844851     ClpObjective * savedObjective = NULL;
     852#ifdef CLP_USEFUL_PRINTOUT
     853     debugInt[0]=numberRows();
     854     debugInt[1]=numberColumns();
     855     debugInt[2]=matrix()->getNumElements();
     856#endif
    845857     if (!objective_ || !matrix_) {
    846858          // totally empty
     
    915927#endif
    916928#endif
     929#ifndef CLPSOLVE_ACTIONS
     930#define CLPSOLVE_ACTIONS 2
     931#endif
     932#if CLPSOLVE_ACTIONS
     933     bool wasAutomatic = (method==ClpSolve::automatic);
     934#endif
    917935     if (presolve != ClpSolve::presolveOff) {
    918936          bool costedSlacks = false;
    919 #ifdef ABC_INHERIT
     937#if CLP_INHERIT_MODE>1
    920938          int numberPasses = 20;
    921939#else
     
    977995            //setSpecialOptions(specialOptions()&~65536);
    978996            // try setting tolerances up
    979 #define CLP_NEW_TOLERANCE 1.0e-7
    980             if (model2->primalTolerance()==1.0e-7&&model2->dualTolerance()==1.0e-7) {
    981               model2->setPrimalTolerance(CLP_NEW_TOLERANCE);
    982               model2->setDualTolerance(CLP_NEW_TOLERANCE);
     997#if CLPSOLVE_ACTIONS
     998            bool changeTolerances=wasAutomatic;
     999#if CLPSOLVE_ACTIONS>1
     1000            changeTolerances=true;
     1001#endif
     1002            if (changeTolerances && model2 != this) {
     1003#define CLP_NEW_TOLERANCE 1.0e-6
     1004              if (model2->primalTolerance()==1.0e-7&&model2->dualTolerance()==1.0e-7) {
     1005                model2->setPrimalTolerance(CLP_NEW_TOLERANCE);
     1006                model2->setDualTolerance(CLP_NEW_TOLERANCE);
     1007              }
    9831008            }
     1009#endif
    9841010#endif
    9851011              model2->eventHandler()->setSimplex(model2);
     
    10101036          }
    10111037     }
     1038#ifdef CLP_USEFUL_PRINTOUT
     1039     debugInt[3]=model2->numberRows();
     1040     debugInt[4]=model2->numberColumns();
     1041     debugInt[5]=model2->matrix()->getNumElements();
     1042     // analyze
     1043     {
     1044       double time1 =CoinCpuTime();
     1045       int numberColumns = model2->numberColumns();
     1046       const double * columnLower = model2->columnLower();
     1047       const double * columnUpper = model2->columnUpper();
     1048       int numberRows = model2->numberRows();
     1049       const double * rowLower = model2->rowLower();
     1050       const double * rowUpper = model2->rowUpper();
     1051       const double * objective = model2->objective();
     1052       CoinPackedMatrix * matrix = model2->matrix();
     1053       CoinBigIndex numberElements = matrix->getNumElements();
     1054       const int * columnLength = matrix->getVectorLengths();
     1055       //const CoinBigIndex * columnStart = matrix->getVectorStarts();
     1056       const double * elementByColumn = matrix->getElements();
     1057       const int * row = matrix->getIndices();
     1058       int * rowCount = new int [numberRows];
     1059       memset(rowCount,0,numberRows*sizeof(int));
     1060       int n=CoinMax (2*numberRows,numberElements);
     1061       n= CoinMax(2*numberColumns,n);
     1062       double * check = new double[n];
     1063       memcpy(check,elementByColumn,numberElements*sizeof(double));
     1064       for (int i=0;i<numberElements;i++) {
     1065         check[i]=fabs(check[i]);
     1066         rowCount[row[i]]++;
     1067       }
     1068       int largestIndex=0;
     1069       for (int i=0;i<numberColumns;i++) {
     1070         largestIndex=CoinMax(largestIndex,columnLength[i]);
     1071       }
     1072       debugInt[12]=largestIndex;
     1073       largestIndex=0;
     1074       for (int i=0;i<numberRows;i++) {
     1075         largestIndex=CoinMax(largestIndex,rowCount[i]);
     1076       }
     1077       n=numberElements;
     1078       delete [] rowCount;
     1079       debugInt[11]=largestIndex;
     1080       std::sort(check,check+n);
     1081       debugDouble[4]=check[0];
     1082       debugDouble[5]=check[n-1];
     1083       int nAtOne=0;
     1084       int nDifferent=0;
     1085       double last=-COIN_DBL_MAX;
     1086       for (int i=0;i<n;i++) {
     1087         if (fabs(last-check[i])>1.0e-12) {
     1088           nDifferent++;
     1089           last=check[i];
     1090         }
     1091         if (check[i]==1.0)
     1092           nAtOne++;
     1093       }
     1094       debugInt[10]=nDifferent;
     1095       debugInt[15]=nAtOne;
     1096       int nInf=0;
     1097       int nZero=0;
     1098       n=0;
     1099       for (int i=0;i<numberRows;i++) {
     1100         double value=fabs(rowLower[i]);
     1101         if (!value)
     1102           nZero++;
     1103         else if (value!=COIN_DBL_MAX)
     1104           check[n++]=value;
     1105         else
     1106           nInf++;
     1107       }
     1108       for (int i=0;i<numberRows;i++) {
     1109         double value=fabs(rowUpper[i]);
     1110         if (!value)
     1111           nZero++;
     1112         else if (value!=COIN_DBL_MAX)
     1113           check[n++]=value;
     1114         else
     1115           nInf++;
     1116       }
     1117       debugInt[16]=nInf;
     1118       debugInt[20]=nZero;
     1119       if (n) {
     1120         std::sort(check,check+n);
     1121         debugDouble[0]=check[0];
     1122         debugDouble[1]=check[n-1];
     1123       } else {
     1124         debugDouble[0]=0.0;
     1125         debugDouble[1]=0.0;
     1126       }
     1127       nAtOne=0;
     1128       nDifferent=0;
     1129       last=-COIN_DBL_MAX;
     1130       for (int i=0;i<n;i++) {
     1131         if (fabs(last-check[i])>1.0e-12) {
     1132           nDifferent++;
     1133           last=check[i];
     1134         }
     1135         if (check[i]==1.0)
     1136           nAtOne++;
     1137       }
     1138       debugInt[8]=nDifferent;
     1139       debugInt[13]=nAtOne;
     1140       nZero=0;
     1141       n=0;
     1142       for (int i=0;i<numberColumns;i++) {
     1143         double value=fabs(objective[i]);
     1144         if (value)
     1145           check[n++]=value;
     1146         else
     1147           nZero++;
     1148       }
     1149       debugInt[21]=nZero;
     1150       if (n) {
     1151         std::sort(check,check+n);
     1152         debugDouble[2]=check[0];
     1153         debugDouble[3]=check[n-1];
     1154       } else {
     1155         debugDouble[2]=0.0;
     1156         debugDouble[3]=0.0;
     1157       }
     1158       nAtOne=0;
     1159       nDifferent=0;
     1160       last=-COIN_DBL_MAX;
     1161       for (int i=0;i<n;i++) {
     1162         if (fabs(last-check[i])>1.0e-12) {
     1163           nDifferent++;
     1164           last=check[i];
     1165         }
     1166         if (check[i]==1.0)
     1167           nAtOne++;
     1168       }
     1169       debugInt[9]=nDifferent;
     1170       debugInt[14]=nAtOne;
     1171       nInf=0;
     1172       nZero=0;
     1173       n=0;
     1174       for (int i=0;i<numberColumns;i++) {
     1175         double value=fabs(columnLower[i]);
     1176         if (!value)
     1177           nZero++;
     1178         else if (value!=COIN_DBL_MAX)
     1179           check[n++]=value;
     1180         else
     1181           nInf++;
     1182       }
     1183       for (int i=0;i<numberColumns;i++) {
     1184         double value=fabs(columnUpper[i]);
     1185         if (!value)
     1186           nZero++;
     1187         else if (value!=COIN_DBL_MAX)
     1188           check[n++]=value;
     1189         else
     1190           nInf++;
     1191       }
     1192       debugInt[17]=nInf;
     1193       double smallestColBound;
     1194       double largestColBound;
     1195       if (n) {
     1196         std::sort(check,check+n);
     1197         smallestColBound=check[0];
     1198         largestColBound=check[n-1];
     1199       } else {
     1200         smallestColBound=0.0;
     1201         largestColBound=0.0;
     1202       }
     1203       nAtOne=0;
     1204       nDifferent=0;
     1205       last=-COIN_DBL_MAX;
     1206       for (int i=0;i<n;i++) {
     1207         if (fabs(last-check[i])>1.0e-12) {
     1208           nDifferent++;
     1209           last=check[i];
     1210         }
     1211         if (check[i]==1.0)
     1212           nAtOne++;
     1213       }
     1214       //debugInt[8]=nDifferent;
     1215       //debugInt[13]=nAtOne;
     1216       printf("BENCHMARK_STATS rhs %d different - %g -> %g (%d at one, %d infinite, %d zero)\n",
     1217              debugInt[8],debugDouble[0],debugDouble[1],debugInt[13],debugInt[16],debugInt[20]);
     1218       printf("BENCHMARK_STATS col %d different - %g -> %g (%d at one, %d infinite, %d zero)\n",
     1219              nDifferent,smallestColBound,largestColBound,nAtOne,nInf,nZero);
     1220       printf("BENCHMARK_STATS els %d different - %g -> %g (%d at one) - longest r,c %d,%d\n",
     1221              debugInt[10],debugDouble[4],debugDouble[5],debugInt[15],
     1222              debugInt[11],debugInt[12]);
     1223       printf("BENCHMARK_STATS obj %d different - %g -> %g (%d at one, %d zero) - time %g\n",
     1224              debugInt[9],debugDouble[2],debugDouble[3],debugInt[14],debugInt[21],
     1225              CoinCpuTime()-time1);
     1226       delete [] check;
     1227     }
     1228#endif
    10121229     if (interrupt)
    10131230          currentModel = model2;
     1231     int saveMoreOptions = moreSpecialOptions_;
    10141232     // For below >0 overrides
    10151233     // 0 means no, -1 means maybe
     
    10201238     int primalStartup = 1;
    10211239     model2->eventHandler()->event(ClpEventHandler::presolveBeforeSolve);
    1022      bool tryItSave = false;
     1240     int tryItSave = 0;
     1241#if CLPSOLVE_ACTIONS
     1242     if (method == ClpSolve::automatic )
     1243       model2->moreSpecialOptions_ |= 8192; // stop switch over
     1244#endif
    10231245     // switch to primal from automatic if just one cost entry
    1024      if (method == ClpSolve::automatic && model2->numberColumns() > 5000 &&
    1025                (specialOptions_ & 1024) != 0) {
     1246     if (method == ClpSolve::automatic && model2->numberColumns() > 5000
     1247#ifndef CLPSOLVE_ACTIONS
     1248         && (specialOptions_ & 1024) != 0
     1249#endif
     1250         ) {
     1251          // look at original model for objective
    10261252          int numberColumns = model2->numberColumns();
    10271253          int numberRows = model2->numberRows();
    1028           const double * obj = model2->objective();
     1254          int numberColumnsOrig = this->numberColumns();
     1255          const double * obj = this->objective();
    10291256          int nNon = 0;
    1030           for (int i = 0; i < numberColumns; i++) {
    1031                if (obj[i])
    1032                     nNon++;
    1033           }
    1034           if (nNon == 1) {
     1257          bool allOnes=true;
     1258          for (int i = 0; i < numberColumnsOrig; i++) {
     1259            if (obj[i]) {
     1260              nNon++;
     1261              if (fabs(obj[i])!=1.0)
     1262                allOnes=false;
     1263            }
     1264          }
     1265          if (nNon <= 1 || allOnes) {
    10351266#ifdef COIN_DEVELOP
    10361267               printf("Forcing primal\n");
    10371268#endif
    10381269               method = ClpSolve::usePrimal;
    1039                tryItSave = numberRows > 200 && numberColumns > 2000 &&
    1040                            (numberColumns > 2 * numberRows || (specialOptions_ & 1024) != 0);
     1270#ifndef CLPSOLVE_ACTIONS
     1271               tryItSave = (numberRows > 200 && numberColumns > 2000 &&
     1272                            (numberColumns > 2 * numberRows || (specialOptions_ & 1024) != 0)) ? 3 : 0;
     1273#else
     1274               if (numberRows > 200 && numberColumns > 2000 &&
     1275                   numberColumns > 2 * numberRows) {
     1276                 tryItSave= 3;
     1277               } else {
     1278                 // If rhs also rubbish then maybe
     1279                 int numberRowsOrig = this->numberRows();
     1280                 const double * rowLower = this->rowLower();
     1281                 const double * rowUpper = this->rowUpper();
     1282                 double last=COIN_DBL_MAX;
     1283                 int nDifferent=0;
     1284                 for (int i=0;i<numberRowsOrig;i++) {
     1285                   double value=fabs(rowLower[i]);
     1286                   if (value&&value!=COIN_DBL_MAX) {
     1287                     if (value!=last) {
     1288                       nDifferent++;
     1289                       last=value;
     1290                     }
     1291                   }
     1292                   value=fabs(rowUpper[i]);
     1293                   if (value&&value!=COIN_DBL_MAX) {
     1294                     if (value!=last) {
     1295                       nDifferent++;
     1296                       last=value;
     1297                     }
     1298                   }
     1299                 }
     1300                 if (nDifferent<2)
     1301                   tryItSave=1;
     1302               }
     1303#endif
    10411304          }
    10421305     }
     
    13361599                         }
    13371600                    }
    1338                     bool tryIt = numberRows > 200 && numberColumns > 2000 &&
    1339                                  (numberColumns > 2 * numberRows || (method != ClpSolve::useDual && (specialOptions_ & 1024) != 0));
     1601                    int tryIt = (numberRows > 200 && numberColumns > 2000 &&
     1602                                 (numberColumns > 2 * numberRows || (method != ClpSolve::useDual && (specialOptions_ & 1024) != 0))) ? 3 : 0;
    13401603                    tryItSave = tryIt;
    13411604                    if (numberRows < 1000 && numberColumns < 3000)
    1342                          tryIt = false;
     1605                         tryIt = 0;
    13431606                    if (notInteger)
    1344                          tryIt = false;
     1607                         tryIt = 0;
    13451608                    if (largest / smallest > 10 || (largest / smallest > 2.0 && largest > 50))
    1346                          tryIt = false;
     1609                         tryIt = 0;
    13471610                    if (tryIt) {
    13481611                         if (largest / smallest > 2.0) {
     
    13821645                    if (!tryIt || nPasses <= 5)
    13831646                         doIdiot = 0;
     1647#if CLPSOLVE_ACTIONS
     1648                    if (doIdiot&&doSprint<0&&wasAutomatic &&
     1649                        20*model2->numberRows()>model2->numberColumns())
     1650                      doSprint=0; // switch off sprint
     1651#endif
    13841652                    if (doSprint) {
    13851653                         method = ClpSolve::usePrimalorSprint;
     
    14441712     }
    14451713     if (method == ClpSolve::useDual) {
     1714#ifdef CLP_USEFUL_PRINTOUT
     1715       debugInt[6]=1;
     1716#endif
    14461717          double * saveLower = NULL;
    14471718          double * saveUpper = NULL;
     
    16601931          timeX = time2;
    16611932     } else if (method == ClpSolve::usePrimal) {
     1933#ifdef CLP_USEFUL_PRINTOUT
     1934       debugInt[6]=2;
     1935#endif
    16621936#ifndef SLIM_CLP
    16631937          if (doIdiot) {
     
    16991973                         increaseSprint = true;
    17001974               }
    1701                bool tryIt = tryItSave;
     1975               int tryIt = tryItSave;
    17021976               if (numberRows < 1000 && numberColumns < 3000)
    1703                     tryIt = false;
     1977                    tryIt = 0;
    17041978               if (tryIt) {
    17051979                    if (increaseSprint) {
     
    17882062                         //info.setStartingWeight(1.0e-1);
    17892063                    }
     2064                    if (tryIt==1)
     2065                      model2->setSpecialOptions(model2->specialOptions()
     2066                                                |8388608);
    17902067               }
    17912068               if (doIdiot > 0) {
     
    18032080                      timeX = time2;
    18042081                    }
     2082#endif
     2083#ifdef CLP_USEFUL_PRINTOUT
     2084                    debugInt[6]=4;
     2085                    debugInt[7]=nPasses;
    18052086#endif
    18062087                    if (nPasses > 70) {
     
    18982179                    info.setStrategy(32 | info.getStrategy());
    18992180                    int saveScalingFlag=model2->scalingFlag();
    1900                     if (objective_->type()==2)
     2181                    bool linearObjective = objective_->type()<2;
     2182                    if (!linearObjective)
    19012183                      model2->scaling(0);
     2184#define CLP_DUAL_IDIOT
     2185#ifdef CLP_DUAL_IDIOT
     2186                    bool doubleIdiot=false;
     2187                    if (nPasses==99&&linearObjective)
     2188                      doubleIdiot=true;
     2189                    if (doubleIdiot) {
     2190                      ClpSimplex * dualModel2 = static_cast<ClpSimplexOther *> (model2)->dualOfModel(1.0,1.0);
     2191                      if (dualModel2) {
     2192                        printf("Dual of model has %d rows and %d columns\n",
     2193                               dualModel2->numberRows(), dualModel2->numberColumns());
     2194                        dualModel2->setOptimizationDirection(1.0);
     2195                        Idiot infoDual(info);
     2196                        infoDual.setModel(dualModel2);
     2197#if 0
     2198                        info.setStrategy(512 | info.getStrategy());
     2199                        // Allow for scaling
     2200                        info.setStrategy(32 | info.getStrategy());
     2201                        info.setStartingWeight(1.0e3);
     2202                        info.setReduceIterations(6);
     2203#endif
     2204                        info.crash(nPasses, model2->messageHandler(),
     2205                                   model2->messagesPointer(),false);
     2206                        infoDual.crash(nPasses, model2->messageHandler(),
     2207                                   model2->messagesPointer(),false);
     2208                        // two copies of solutions
     2209                        ClpSimplex temp(*model2);
     2210#if 0
     2211                        static_cast<ClpSimplexOther *> (&temp)->restoreFromDual(dualModel2);
     2212#else
     2213                        // move duals and just copy primal
     2214                        memcpy(temp.dualRowSolution(),dualModel2->primalColumnSolution(),
     2215                               numberRows*sizeof(double));
     2216                        memcpy(temp.primalColumnSolution(),model2->primalColumnSolution(),
     2217                               numberColumns*sizeof(double));
     2218#endif
     2219                        delete dualModel2;
     2220                        int numberRows=model2->numberRows();
     2221                        int numberColumns=model2->numberColumns();
     2222                        ClpSimplex * tempModel[2];
     2223                        tempModel[0]=model2;
     2224                        tempModel[1]=&temp;
     2225                        const double * primalColumn[2];
     2226                        const double * dualRow[2];
     2227                        double * dualColumn[2];
     2228                        double * primalRow[2];
     2229                        for (int i=0;i<2;i++) {
     2230                          primalColumn[i]=tempModel[i]->primalColumnSolution();
     2231                          dualRow[i]=tempModel[i]->dualRowSolution();
     2232                          dualColumn[i]=tempModel[i]->dualColumnSolution();
     2233                          primalRow[i]=tempModel[i]->primalRowSolution();
     2234                          memcpy(dualColumn[i],model2->objective(),
     2235                                 numberColumns*sizeof(double));
     2236                          memset(primalRow[i], 0, numberRows * sizeof(double));
     2237                          tempModel[i]->clpMatrix()->transposeTimes(-1.0,
     2238                                                              dualRow[i],
     2239                                                              dualColumn[i]);
     2240                          tempModel[i]->clpMatrix()->times(1.0,
     2241                                                     primalColumn[i],
     2242                                                     primalRow[i]);
     2243                          tempModel[i]->checkSolutionInternal();
     2244                          printf("model %d - dual inf %g primal inf %g\n",
     2245                                 i,tempModel[i]->sumDualInfeasibilities(),
     2246                                 tempModel[i]->sumPrimalInfeasibilities());
     2247                        }
     2248                        printf("What now\n");
     2249                      } else {
     2250                        doubleIdiot=false;
     2251                      }
     2252                    }
     2253                    if (!doubleIdiot)
     2254                      info.crash(nPasses, model2->messageHandler(), model2->messagesPointer(),
     2255                                (objective_->type() <2));
     2256#else
    19022257                    info.crash(nPasses, model2->messageHandler(), model2->messagesPointer(),(objective_->type() <2));
     2258#endif
    19032259                      model2->scaling(saveScalingFlag);
    19042260#endif
     
    22532609          // but if big use to get ratio
    22542610          double ratio = 3;
     2611#ifdef CLP_USEFUL_PRINTOUT
     2612          debugInt[6]=3;
     2613          debugInt[7]=maxSprintPass;
     2614#endif
    22552615          if (maxSprintPass > 1000) {
    22562616               ratio = static_cast<double> (maxSprintPass) * 0.0001;
     
    32293589          int rcode=eventHandler()->event(ClpEventHandler::presolveAfterFirstSolve);
    32303590          if (finalStatus != 3 && rcode < 0 && (finalStatus || oldStatus == -1)) {
     3591               double sumPrimal=sumPrimalInfeasibilities_;
     3592               double sumDual=sumDualInfeasibilities_;
     3593               // ignore some parts of solution
     3594               if (finalStatus == 1) {
     3595                 // infeasible
     3596                 sumDual=0.0;
     3597               } else if (finalStatus == 2) {
     3598                 sumPrimal=0.0;
     3599               }
    32313600               int savePerturbation = perturbation();
    3232                if (!finalStatus || (moreSpecialOptions_ & 2) == 0 ||
    3233                    fabs(sumDualInfeasibilities_)+
    3234                    fabs(sumPrimalInfeasibilities_)<1.0e-3) {
     3601               if (savePerturbation==50)
     3602                 setPerturbation(51); // small
     3603               if (!finalStatus || finalStatus == 2 ||(moreSpecialOptions_ & 2) == 0 ||
     3604                   fabs(sumDual)+
     3605                   fabs(sumPrimal)<1.0e-3) {
    32353606                    if (finalStatus == 2) {
     3607                      if (sumDual>1.0e-4) {
    32363608                         // unbounded - get feasible first
    32373609                         double save = optimizationDirection_;
     
    32393611                         primal(1);
    32403612                         optimizationDirection_ = save;
    3241                          primal(1);
     3613                      }
     3614                      primal(1);
    32423615                    } else if (finalStatus == 1) {
    32433616                         dual();
    32443617                    } else {
    32453618                        if ((moreSpecialOptions_&65536)==0) {
    3246                           if (numberRows_<10000)
     3619                          if (numberRows_<10000||true)
    32473620                            setPerturbation(100); // probably better to perturb after n its
    32483621                          else if (savePerturbation<100)
     
    33353708     eventHandler()->event(ClpEventHandler::presolveEnd);
    33363709     delete pinfo;
     3710     moreSpecialOptions_= saveMoreOptions;
    33373711     return finalStatus;
    33383712}
  • trunk/Clp/src/ClpSolve.hpp

    r2025 r2074  
    204204     /// Whether we want to kill small substitutions
    205205     inline bool doKillSmall() const {
    206           return (independentOptions_[1] & 1024) == 0;
     206          return (independentOptions_[1] & 8192) == 0;
    207207     }
    208208     inline void setDoKillSmall(bool doKill) {
    209           if (doKill) independentOptions_[1]  &= ~1024;
    210           else independentOptions_[1] |= 1024;
     209          if (doKill) independentOptions_[1]  &= ~8192;
     210          else independentOptions_[1] |= 8192;
    211211     }
    212212     /// Set whole group
  • trunk/Clp/src/ClpSolver.cpp

    r2070 r2074  
    7373#include "dmalloc.h"
    7474#endif
     75#ifdef CLP_USEFUL_PRINTOUT
     76static double startElapsed=0.0;
     77static double startCpu=0.0;
     78static std::string mpsFile="";
     79extern double debugDouble[10];
     80extern int debugInt[24];
     81#endif
    7582#if defined(COIN_HAS_WSMP) || defined(COIN_HAS_AMD) || defined(COIN_HAS_CHOLMOD) || defined(TAUCS_BARRIER) || defined(COIN_HAS_MUMPS)
    7683#define FOREIGN_BARRIER
     
    97104          return;
    98105     }
     106  void openblas_set_num_threads(int num_threads);
    99107}
    100108
     
    161169          models->setPerturbation(50);
    162170          models->messageHandler()->setPrefix(false);
    163 #ifdef ABC_INHERIT
     171#if CLP_INHERIT_MODE>1
    164172          models->setDualTolerance(1.0e-6);
    165173          models->setPrimalTolerance(1.0e-6);
     
    175183          // Set up all non-standard stuff
    176184          //int numberModels=1;
     185#ifdef CLP_USEFUL_PRINTOUT
     186          startElapsed=CoinGetTimeOfDay();
     187          startCpu=CoinCpuTime();
     188          memset(debugInt,0,sizeof(debugInt));
     189          memset(debugDouble,0,sizeof(debugDouble));
     190#endif
    177191#if 0
    178192#ifndef ABC_INHERIT
     
    196210          int doSprint = -1;
    197211          // set reasonable defaults
    198 #ifdef ABC_INHERIT
     212#if CLP_INHERIT_MODE>1
    199213#define DEFAULT_PRESOLVE_PASSES 20
    200214#else
     
    810824                         case CBC_PARAM_ACTION_BAB:
    811825                              if (goodModels[iModel]) {
     826                                //openblas_set_num_threads(4);
    812827                                if (type==CLP_PARAM_ACTION_EITHERSIMPLEX||
    813                                     type==CBC_PARAM_ACTION_BAB)
     828                                    type==CBC_PARAM_ACTION_BAB)
    814829                                  models[iModel].setMoreSpecialOptions(16384|
    815830                                                                       models[iModel].moreSpecialOptions());
     
    12201235                                             currentModel = NULL;
    12211236                                        }
     1237#ifndef COIN_HAS_ASL
     1238                                        CoinMessageHandler * generalMessageHandler = models->messageHandler();
     1239                                        generalMessageHandler->setPrefix(false);
     1240                                        CoinMessages generalMessages = models->messages();
     1241                                        char generalPrint[100];
     1242#endif
     1243                                        sprintf(generalPrint, "After translating dual back to primal - objective value is %g",
     1244                                                thisModel->objectiveValue());
     1245                                        generalMessageHandler->message(CLP_GENERAL, generalMessages)
     1246                                        << generalPrint
     1247                                        << CoinMessageEol;
    12221248                                        // switch off (user can switch back on)
    12231249                                        parameters[whichParam(CLP_PARAM_INT_DUALIZE,
     
    14131439                              if (canOpen) {
    14141440                                   int status;
     1441#ifdef CLP_USEFUL_PRINTOUT
     1442                                   mpsFile=fileName;
     1443#endif
    14151444                                   if (!gmpl)
    14161445                                        status = models[iModel].readMps(fileName.c_str(),
     
    28752904     dmalloc_log_unfreed();
    28762905     dmalloc_shutdown();
     2906#endif
     2907#ifdef CLP_USEFUL_PRINTOUT
     2908     printf("BENCHMARK %s took %g cpu seconds (%g elapsed) - %d row, %d columns, %d elements - reduced to %d, %d, %d\n",
     2909              mpsFile.c_str(),CoinCpuTime()-startCpu,CoinGetTimeOfDay()-startElapsed,
     2910              debugInt[0],debugInt[1],debugInt[2],
     2911              debugInt[3],debugInt[4],debugInt[5]);
     2912            const char * method[]={"","Dual","Primal","Sprint","Idiot"};
     2913            printf("BENCHMARK using method %s (%d)\n",
     2914              method[debugInt[6]],debugInt[7]);
    28772915#endif
    28782916     return 0;
  • trunk/Clp/src/CoinAbcBaseFactorization1.cpp

    r2015 r2074  
    24492449#else
    24502450        status=2;
    2451 #if 0 //ABC_NORMAL_DEBUG>1
     2451#if 1 //0 //ABC_NORMAL_DEBUG>1
    24522452        std::cout<<"      Went dense at "<<numberRowsLeft_<<" rows "<<
    24532453          totalElements_<<" "<<full<<" "<<leftElements<<std::endl;
  • trunk/Clp/src/CoinAbcCommon.hpp

    r1910 r2074  
    4141#ifndef FAKE_CILK
    4242#include <cilk/cilk.h>
     43//gcc4.9 main branch has not got cilk_for!!!!!!!!!!!
     44#ifdef GCC_4_9
     45#undef cilk_for
     46#define cilk_for for
     47#endif
    4348#else
    4449#define cilk_for for
     
    191196enum CBLAS_TRANSPOSE {CblasNoTrans=111, CblasTrans=112, CblasConjTrans=113,
    192197                      AtlasConj=114};
    193 #define CLAPACK
     198//#define CLAPACK
    194199// using simple lapack interface
    195200extern "C"
  • trunk/Clp/src/IdiSolve.cpp

    r2030 r2074  
    338338          djSave[i] = 1.0e30;
    339339     }
     340#ifndef OSI_IDIOT
     341     int numberColumns = model_->numberColumns();
     342     for (int i=0;i<numberColumns;i++) {
     343       if (model_->getColumnStatus(i)!=ClpSimplex::isFixed)
     344         statusSave[i] = 0;
     345       else
     346         statusSave[i] = 2;
     347     }
     348     memset(statusSave+numberColumns,0,ncols-numberColumns);
     349     for (int i=0;i<numberColumns;i++) {
     350       if (model_->getColumnStatus(i)==ClpSimplex::isFixed) {
     351         assert (colsol[i]<lower[i]+tolerance||
     352                 colsol[i]>upper[i]-tolerance);
     353       }
     354     }
     355#else
    340356     for (i = 0; i < ncols; i++) {
    341357          if (upper[i] - lower[i]) {
     
    345361          }
    346362     }
     363#endif
    347364     // for two pass method
    348365     int start[2];
  • trunk/Clp/src/Idiot.cpp

    r2070 r2074  
    1313#include "CoinTime.hpp"
    1414#include "CoinSort.hpp"
     15#include "CoinFactorization.hpp"
    1516#include "CoinMessageHandler.hpp"
    1617#include "CoinHelperFunctions.hpp"
     
    475476     pi = new double[nrows];
    476477     dj = new double[ncols];
     478#ifndef OSI_IDIOT
     479     bool fixAfterSome = false; //(model_->specialOptions()&8388608)!=0;
     480     int exitAfter = 1000000; //(model_->specialOptions()&8388608)!=0 ? 50 : 1000000;
     481     {
     482        int numberColumns = model_->numberColumns();
     483        for (int i=0;i<numberColumns;i++) {
     484          if (upper[i]==lower[i])
     485            model_->setColumnStatus(i,ClpSimplex::isFixed);
     486        }
     487     }
     488#endif
    477489     delete [] whenUsed_;
    478490     bool oddSlacks = false;
     
    787799               break;
    788800          iteration++;
     801          if (iteration>=exitAfter)
     802            break;
    789803          checkIteration++;
    790804          if (lastResult.infeas <= smallInfeas && lastResult.objval < bestFeasible) {
     
    835849#endif
    836850          }
     851#ifndef OSI_IDIOT
     852          if (fixAfterSome) {
     853            if (result.infeas<0.01*nrows&&iteration>10&&(3*n>2*nrows||4*n>2*ncols)) {
     854                // flag
     855                int numberColumns = model_->numberColumns();
     856                printf("Flagging satisfied\n");
     857                fixAfterSome=false;
     858                for (int i=0;i<numberColumns;i++) {
     859                  if (colsol[i]>upper[i]-1.0e-7||
     860                      colsol[i]<lower[i]+1.0e-7)
     861                    model_->setColumnStatus(i,ClpSimplex::isFixed);
     862                }
     863              }
     864              }
     865#endif
    837866          if (iteration > 50 && n == numberAway ) {
    838867            if((result.infeas < 1.0e-4 && majorIterations_<200)||result.infeas<1.0e-8) {
     
    13061335     }
    13071336#endif
     1337#define FEB_TRY
    13081338#ifdef FEB_TRY
    13091339     int savePerturbation = model_->perturbation();
    13101340     int saveOptions = model_->specialOptions();
    13111341     model_->setSpecialOptions(saveOptions | 8192);
    1312      if (savePerturbation_ == 50)
    1313           model_->setPerturbation(56);
     1342     //if (savePerturbation_ == 50)
     1343     //   model_->setPerturbation(56);
    13141344#endif
    13151345     model_->createStatus();
     
    16091639     } else {
    16101640          /* still try and put singletons rather than artificials in basis */
    1611           int ninbas = 0;
    16121641          for (i = 0; i < nrows; i++) {
    1613                model_->setRowStatus(i, ClpSimplex::basic);
    1614           }
    1615           for (i = 0; i < ncols; i++) {
    1616                if (columnLength[i] == 1 && upper[i] > lower[i] + 1.0e-5) {
     1642              model_->setRowStatus(i, ClpSimplex::basic);
     1643          }
     1644          int ninbas = 0;
     1645          for (i = 0; i < ncols; i++) {
     1646              if (columnLength[i] == 1 && upper[i] > lower[i] + 1.0e-5) {
    16171647                    CoinBigIndex j = columnStart[i];
    16181648                    double value = element[j];
     
    16761706                    }
    16771707               }
    1678           }
    1679           /*printf("%d in basis\n",ninbas);*/
     1708            }
     1709            /*printf("%d in basis\n",ninbas);*/
    16801710     }
    16811711     bool wantVector = false;
     
    16851715          wantVector = clpMatrixO->wantsSpecialColumnCopy();
    16861716     }
     1717     if ((model_->specialOptions()&8388608)!=0) { // temp test
     1718       justValuesPass=true;
     1719       assert (addAll<3);
     1720       assert (presolve);
     1721       allowInfeasible=true;
     1722     }
     1723     double * saveBounds=NULL;
    16871724     if (addAll < 3) {
    16881725          ClpPresolve pinfo;
     
    16901727               if (allowInfeasible) {
    16911728                    // fix up so will be feasible
    1692                     double * rhs = new double[nrows];
     1729                    const double * dual = model_->dualRowSolution();
     1730                    double * rhs = new double[nrows];
     1731                    double * saveBounds=new double[2*ncols];
     1732                    memcpy(saveBounds,lower,ncols*sizeof(double));
     1733                    memcpy(saveBounds+ncols,upper,ncols*sizeof(double));
     1734                    for (i = 0; i < nrows; i++)
     1735                      rhs[i]=fabs(dual[i]);
     1736                    std::sort(rhs,rhs+nrows);
     1737                    int nSmall=nrows;
     1738                    int nMedium=nrows;
     1739                    double small=CoinMax(1.0e-4,1.0e-5*rhs[nrows-1]);
     1740                    double medium=CoinMax(1.0e-2,1.0e-3*rhs[nrows-1]);
     1741                    char * marked = new char [nrows];
     1742                    double * rowupper = model_->rowUpper();
     1743                    double * rowlower = model_->rowLower();
     1744                    // if tiny then drop row??
     1745                    for (i = 0; i < nrows; i++) {
     1746                      if (rhs[i]>=small) {
     1747                        nSmall=i-1;
     1748                        break;
     1749                      }
     1750                    }
     1751                    for (; i < nrows; i++) {
     1752                      if (rhs[i]>=medium) {
     1753                        nMedium=i-1;
     1754                        break;
     1755                      }
     1756                    }
     1757                    printf("%d < %g, %d < %g, %d >= %g\n",
     1758                           nSmall,small,nMedium-nSmall,medium,nrows-nMedium,medium);
     1759                    for (i = 0; i < nrows; i++) {
     1760                      if (fabs(dual[i])<medium) {
     1761                        marked[i]=0;
     1762                      } else {
     1763                        marked[i]=1;
     1764                      }
     1765                    }
    16931766                    memset(rhs, 0, nrows * sizeof(double));
    16941767                    model_->clpMatrix()->times(1.0, colsol, rhs);
    1695                     double * rowupper = model_->rowUpper();
    1696                     double * rowlower = model_->rowLower();
    16971768                    saveRowUpper = CoinCopyOfArray(rowupper, nrows);
    16981769                    saveRowLower = CoinCopyOfArray(rowlower, nrows);
    16991770                    double sum = 0.0;
     1771                    int nFixedRows=0;
    17001772                    for (i = 0; i < nrows; i++) {
     1773                         if (marked[i]&&rowupper[i]>rowlower[i]) {
     1774                           double distance=CoinMin(rowupper[i]-rhs[i],rhs[i]-rowlower[i]);
     1775                           // look at distance and sign of dual
     1776                           if (rhs[i] > rowupper[i]-1.0e-6) {
     1777                             rhs[i]=rowupper[i];
     1778                             rowlower[i]=rowupper[i];
     1779                             nFixedRows++;
     1780                           } else if (rhs[i] < rowlower[i]+1.0e-6) {
     1781                             rhs[i]=rowlower[i];
     1782                             rowupper[i]=rowlower[i];
     1783                             nFixedRows++;
     1784                           }
     1785                         }
    17011786                         if (rhs[i] > rowupper[i]) {
    17021787                              sum += rhs[i] - rowupper[i];
    17031788                              rowupper[i] = rhs[i];
     1789                              // maybe make equality
    17041790                         }
    17051791                         if (rhs[i] < rowlower[i]) {
    17061792                              sum += rowlower[i] - rhs[i];
    17071793                              rowlower[i] = rhs[i];
    1708                          }
    1709                     }
    1710                     COIN_DETAIL_PRINT(printf("sum of infeasibilities %g\n", sum));
     1794                              // maybe make equality
     1795                         }
     1796                    }
     1797                    int nFixed=0;
     1798                    for (i = 0; i < ncols; i++) {
     1799                      if (colsol[i]<lower[i]+1.0e-7) {
     1800                        upper[i]=lower[i];
     1801                        nFixed++;
     1802                      } else if (colsol[i]>upper[i]-1.0e-7) {
     1803                        lower[i]=upper[i];
     1804                        nFixed++;
     1805                      }
     1806                    }
     1807                    printf("sum of infeasibilities %g - %d fixed rows, %d fixed columns\n",
     1808                           sum,nFixedRows,nFixed);
     1809                    //COIN_DETAIL_PRINT(printf("sum of infeasibilities %g\n", sum));
    17111810                    delete [] rhs;
     1811                    delete [] marked;
    17121812               }
    17131813               saveModel = model_;
    17141814               pinfo.setPresolveActions(pinfo.presolveActions() | 16384);
    17151815               model_ = pinfo.presolvedModel(*model_, 1.0e-8, false, 5);
     1816               if (saveBounds) {
     1817                 memcpy(saveModel->columnLower(),saveBounds,ncols*sizeof(double));
     1818                 memcpy(saveModel->columnUpper(),saveBounds+ncols,ncols*sizeof(double));
     1819                 delete [] saveBounds;
     1820               }
     1821               if (model_&&(saveModel->specialOptions()&8388608)!=0) {
     1822                 int nrows = model_->getNumRows();
     1823                 int ncols = model_->getNumCols();
     1824                 double * lower = model_->columnLower();
     1825                 double * upper = model_->columnUpper();
     1826                 const double * rowlower = model_->getRowLower();
     1827                 const double * rowupper = model_->getRowUpper();
     1828                 double * rowsol = model_->primalRowSolution();
     1829                 double * colsol = model_->primalColumnSolution();;
     1830                 int ninbas = 0;
     1831                 int * which = new int[2*ncols+nrows];
     1832                 double * dj = model_->dualColumnSolution();
     1833                 for (int i=0;i<ncols;i++) {
     1834                   dj[i]=-CoinMin(upper[i]-colsol[i],colsol[i]-lower[i]);
     1835                   which[i]=i;
     1836                 }
     1837                 CoinSort_2(dj,dj+ncols,which);
     1838                 ninbas=CoinMin(ncols,nrows);
     1839                 int * columnIsBasic=which+ncols;
     1840                 int * rowIsBasic=columnIsBasic+ncols;
     1841                 for (int i=0;i<nrows+ncols;i++)
     1842                   columnIsBasic[i]=-1;
     1843                 for (int i=0;i<ninbas;i++) {
     1844                   int iColumn=which[i];
     1845                   columnIsBasic[iColumn]=i;
     1846                 }
     1847                 // factorize
     1848                 CoinFactorization factor;
     1849                 factor.pivotTolerance(0.1);
     1850                 factor.setDenseThreshold(0);
     1851                 int status=-1;
     1852                 // If initial is too dense - then all-slack may be better
     1853                 double areaFactor=4.0;
     1854                 const CoinPackedMatrix * matrix = model_->matrix();
     1855                 while (status) {
     1856                   status=factor.factorize(*matrix,rowIsBasic,columnIsBasic,areaFactor);
     1857                   if (status==-99) {
     1858                     // put all slacks in
     1859                     for (int i=0;i<nrows;i++)
     1860                       rowIsBasic[i]=i;
     1861                     for (int i=0;i<ncols;i++)
     1862                      columnIsBasic[i]=-1;
     1863                     break;
     1864                   } else if (status==-1) {
     1865                     factor.pivotTolerance(0.99);
     1866                     // put all slacks in
     1867                     for (int i=0;i<nrows;i++)
     1868                       rowIsBasic[i]=i;
     1869                     for (int i=0;i<ncols;i++) {
     1870                       int iRow=columnIsBasic[i];
     1871                       if (iRow>=0)
     1872                         rowIsBasic[iRow]=-1; // out
     1873                     }
     1874                   }
     1875                 }
     1876                 delete [] which;
     1877                 for (int i = 0; i < nrows; i++) {
     1878                   if (rowIsBasic[i]>=0) {
     1879                     model_->setRowStatus(i, ClpSimplex::basic);
     1880                   } else if (rowlower[i]==rowupper[i]) {
     1881                     model_->setRowStatus(i, ClpSimplex::isFixed);
     1882                   } else if (rowsol[i]-rowlower[i]<rowupper[i]-rowsol[i]) {
     1883                     model_->setRowStatus(i, ClpSimplex::atLowerBound);
     1884                   } else {
     1885                     model_->setRowStatus(i, ClpSimplex::atUpperBound);
     1886                   }
     1887                 }
     1888                 for (int i=0;i<ncols;i++) {
     1889                   if (colsol[i]>upper[i]-1.0e-7||
     1890                       colsol[i]<lower[i]+1.0e-7) {
     1891                     model_->setColumnStatus(i,ClpSimplex::isFixed);
     1892                   } else if (columnIsBasic[i]>=0) {
     1893                     model_->setColumnStatus(i,ClpSimplex::basic);
     1894                   } else {
     1895                     model_->setColumnStatus(i,ClpSimplex::superBasic);
     1896                   }
     1897                 }
     1898               }
    17161899          }
    17171900          if (model_) {
     
    17441927               }
    17451928               if (presolve) {
    1746                  model_->primal(1);
     1929                 if (!justValuesPass)
     1930                   model_->primal(1);
    17471931                    pinfo.postsolve(true);
    17481932                    delete model_;
    17491933                    model_ = saveModel;
    17501934                    saveModel = NULL;
    1751                }
     1935               } 
    17521936          } else {
    17531937               // not feasible
  • trunk/Clp/src/Idiot.hpp

    r1665 r2074  
    205205          dropEnoughWeighted_ = value;
    206206     }
     207     /// Set model
     208     inline void setModel(OsiSolverInterface * model) {
     209       model_ = model;
     210     };
    207211     //@}
    208212
Note: See TracChangeset for help on using the changeset viewer.