Ignore:
Timestamp:
Jun 26, 2007 3:00:48 AM (13 years ago)
Author:
forrest
Message:

moving branches/devel to trunk

Location:
trunk
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • trunk

    • Property svn:externals
      •  

        old new  
        11MSVisualStudio   https://projects.coin-or.org/svn/MSVisualStudio/trunk/ExternalsDirs/Clp
        22BuildTools    https://projects.coin-or.org/svn/BuildTools/trunk
        3 Data/Netlib   https://projects.coin-or.org/svn/Data/trunk/Netlib
        4 Data/Sample   https://projects.coin-or.org/svn/Data/trunk/Sample
        5 CoinUtils     https://projects.coin-or.org/svn/CoinUtils/stable/2.0/CoinUtils
         3ThirdParty/Blas https://projects.coin-or.org/svn/BuildTools/ThirdParty/Blas/stable/1.0
         4ThirdParty/Lapack https://projects.coin-or.org/svn/BuildTools/ThirdParty/Lapack/stable/1.0
         5Data/Netlib   https://projects.coin-or.org/svn/Data/stable/1.0/Netlib
         6Data/Sample   https://projects.coin-or.org/svn/Data/stable/1.0/Sample
         7CoinUtils     https://projects.coin-or.org/svn/CoinUtils/trunk/CoinUtils
  • trunk/Clp/src/Idiot.cpp

    r1015 r1034  
    1010#include "Idiot.hpp"
    1111#include "CoinTime.hpp"
     12#include "CoinSort.hpp"
    1213#include "CoinMessageHandler.hpp"
    1314#include "CoinHelperFunctions.hpp"
     
    6061  }
    6162}
     63// Deals with whenUsed and slacks
     64int
     65Idiot::cleanIteration(int iteration, int ordinaryStart, int ordinaryEnd,
     66                      double * colsol, const double * lower, const double * upper,
     67                      const double * rowLower, const double * rowUpper,
     68                      const double * cost, const double * element, double fixTolerance,
     69                      double & objValue, double & infValue)
     70{
     71  int n=0;
     72  if ((strategy_&16384)==0) {
     73    for (int i=ordinaryStart;i<ordinaryEnd;i++) {
     74      if (colsol[i]>lower[i]+fixTolerance) {
     75        if (colsol[i]<upper[i]-fixTolerance) {
     76          n++;
     77        } else {
     78          colsol[i]=upper[i];
     79        }
     80        whenUsed_[i]=iteration;
     81      } else {
     82        colsol[i]=lower[i];
     83      }
     84    }
     85    return n;
     86  } else {
     87    printf("entering inf %g, obj %g\n",infValue,objValue);
     88    int nrows=model_->getNumRows();
     89    int ncols=model_->getNumCols();
     90    int * posSlack = whenUsed_+ncols;
     91    int * negSlack = posSlack+nrows;
     92    int * nextSlack = negSlack + nrows;
     93    double * rowsol = (double *) (nextSlack+ncols);
     94    memset(rowsol,0,nrows*sizeof(double));
     95#ifdef OSI_IDIOT
     96    const CoinPackedMatrix * matrix = model_->getMatrixByCol();
     97#else
     98    ClpMatrixBase * matrix = model_->clpMatrix();
     99#endif
     100    const int * row = matrix->getIndices();
     101    const CoinBigIndex * columnStart = matrix->getVectorStarts();
     102    const int * columnLength = matrix->getVectorLengths();
     103    //const double * element = matrix->getElements();
     104    int i;
     105    objValue=0.0;
     106    infValue=0.0;
     107    for ( i=0;i<ncols;i++) {
     108      if (nextSlack[i]==-1) {
     109        // not a slack
     110        if (colsol[i]>lower[i]+fixTolerance) {
     111          if (colsol[i]<upper[i]-fixTolerance) {
     112            n++;
     113            whenUsed_[i]=iteration;
     114          } else {
     115            colsol[i]=upper[i];
     116          }
     117          whenUsed_[i]=iteration;
     118        } else {
     119          colsol[i]=lower[i];
     120        }
     121        double value = colsol[i];
     122        if (value) {
     123          objValue += cost[i]*value;
     124          CoinBigIndex j;
     125          for (j=columnStart[i];j<columnStart[i]+columnLength[i];j++) {
     126            int iRow = row[j];
     127            rowsol[iRow] += value*element[j];
     128          }
     129        }
     130      }
     131    }
     132    // temp fix for infinite lbs - just limit to -1000
     133    for (i=0;i<nrows;i++) {
     134      double rowSave=rowsol[i];
     135      int iCol;
     136      iCol =posSlack[i];
     137      if (iCol>=0) {
     138        // slide all slack down
     139        double rowValue=rowsol[i];
     140        CoinBigIndex j=columnStart[iCol];
     141        double lowerValue = CoinMax(CoinMin(colsol[iCol],0.0)-1000.0,lower[iCol]);
     142        rowSave += (colsol[iCol]-lowerValue)*element[j];
     143        colsol[iCol]=lowerValue;
     144        while (nextSlack[iCol]>=0) {
     145          iCol = nextSlack[iCol];
     146          double lowerValue = CoinMax(CoinMin(colsol[iCol],0.0)-1000.0,lower[iCol]);
     147          j=columnStart[iCol];
     148          rowSave += (colsol[iCol]-lowerValue)*element[j];
     149          colsol[iCol]=lowerValue;
     150        }
     151        iCol =posSlack[i];
     152        while (rowValue<rowLower[i]&&iCol>=0) {
     153          // want to increase
     154          double distance = rowLower[i]-rowValue;
     155          double value = element[columnStart[iCol]];
     156          double thisCost = cost[iCol];
     157          if (distance<=value*(upper[iCol]-colsol[iCol])) {
     158            // can get there
     159            double movement = distance/value;
     160            objValue += movement*thisCost;
     161            rowValue = rowLower[i];
     162            colsol[iCol] += movement;
     163          } else {
     164            // can't get there
     165            double movement = upper[iCol]-colsol[iCol];
     166            objValue += movement*thisCost;
     167            rowValue += movement*value;
     168            colsol[iCol] = upper[iCol];
     169            iCol = nextSlack[iCol];
     170          }
     171        }
     172        if (iCol>=0) {
     173          // may want to carry on - because of cost?
     174          while (cost[iCol]<0&&rowValue<rowUpper[i]) {
     175            // want to increase
     176            double distance = rowUpper[i]-rowValue;
     177            double value = element[columnStart[iCol]];
     178            double thisCost = cost[iCol];
     179            if (distance<=value*(upper[iCol]-colsol[iCol])) {
     180              // can get there
     181              double movement = distance/value;
     182              objValue += movement*thisCost;
     183              rowValue = rowUpper[i];
     184              colsol[iCol] += movement;
     185              iCol=-1;
     186            } else {
     187              // can't get there
     188              double movement = upper[iCol]-colsol[iCol];
     189              objValue += movement*thisCost;
     190              rowValue += movement*value;
     191              colsol[iCol] = upper[iCol];
     192              iCol = nextSlack[iCol];
     193            }
     194          }
     195          if (iCol>=0&&colsol[iCol]>lower[iCol]+fixTolerance&&
     196              colsol[iCol]<upper[iCol]-fixTolerance) {
     197            whenUsed_[i]=iteration;
     198            n++;
     199          }
     200        }
     201        rowsol[i]=rowValue;
     202      }
     203      iCol =negSlack[i];
     204      if (iCol>=0) {
     205        // slide all slack down
     206        double rowValue=rowsol[i];
     207        CoinBigIndex j=columnStart[iCol];
     208        rowSave += (colsol[iCol]-lower[iCol])*element[j];
     209        colsol[iCol]=lower[iCol];
     210        assert (lower[iCol]>-1.0e20);
     211        while (nextSlack[iCol]>=0) {
     212          iCol = nextSlack[iCol];
     213          j=columnStart[iCol];
     214          rowSave += (colsol[iCol]-lower[iCol])*element[j];
     215          colsol[iCol]=lower[iCol];
     216        }
     217        iCol =negSlack[i];
     218        while (rowValue>rowUpper[i]&&iCol>=0) {
     219          // want to increase
     220          double distance = -(rowUpper[i]-rowValue);
     221          double value = -element[columnStart[iCol]];
     222          double thisCost = cost[iCol];
     223          if (distance<=value*(upper[iCol]-lower[iCol])) {
     224            // can get there
     225            double movement = distance/value;
     226            objValue += movement*thisCost;
     227            rowValue = rowUpper[i];
     228            colsol[iCol] += movement;
     229          } else {
     230            // can't get there
     231            double movement = upper[iCol]-lower[iCol];
     232            objValue += movement*thisCost;
     233            rowValue -= movement*value;
     234            colsol[iCol] = upper[iCol];
     235            iCol = nextSlack[iCol];
     236          }
     237        }
     238        if (iCol>=0) {
     239          // may want to carry on - because of cost?
     240          while (cost[iCol]<0&&rowValue>rowLower[i]) {
     241            // want to increase
     242            double distance = -(rowLower[i]-rowValue);
     243            double value = -element[columnStart[iCol]];
     244            double thisCost = cost[iCol];
     245            if (distance<=value*(upper[iCol]-colsol[iCol])) {
     246              // can get there
     247              double movement = distance/value;
     248              objValue += movement*thisCost;
     249              rowValue = rowLower[i];
     250              colsol[iCol] += movement;
     251              iCol=-1;
     252            } else {
     253              // can't get there
     254              double movement = upper[iCol]-colsol[iCol];
     255              objValue += movement*thisCost;
     256              rowValue -= movement*value;
     257              colsol[iCol] = upper[iCol];
     258              iCol = nextSlack[iCol];
     259            }
     260          }
     261          if (iCol>=0&&colsol[iCol]>lower[iCol]+fixTolerance&&
     262              colsol[iCol]<upper[iCol]-fixTolerance) {
     263            whenUsed_[i]=iteration;
     264            n++;
     265          }
     266        }
     267        rowsol[i]=rowValue;
     268      }
     269      infValue += CoinMax(CoinMax(0.0,rowLower[i]-rowsol[i]),rowsol[i]-rowUpper[i]);
     270      // just change
     271      rowsol[i] -= rowSave;
     272    }
     273    return n;
     274  }
     275}
    62276
    63 /* returns -1 or start of costed slacks */
     277/* returns -1 if none or start of costed slacks or -2 if
     278   there are costed slacks but it is messy */
    64279 static int countCostedSlacks(OsiSolverInterface * model)
    65280{
     
    121336  }
    122337  sum /= (double) (nnzero+1);
    123   maxIts_=2;
     338  if (maxIts_==5)
     339    maxIts_=2;
    124340  if (numberPass<=0)
    125341    // Cast to double to avoid VACPP complaining
     
    130346  if (mu_==1e-4)
    131347    mu_= CoinMax(1.0e-3,sum*1.0e-5);
    132   if (!lightWeight_) {
    133     maxIts2_=105;
    134   } else if (lightWeight_==1) {
    135     mu_ *= 1000.0;
    136     maxIts2_=23;
    137   } else if (lightWeight_==2) {
    138     maxIts2_=11;
    139   } else {
    140     maxIts2_=23;
     348  if (maxIts2_==100) {
     349    if (!lightWeight_) {
     350      maxIts2_=105;
     351    } else if (lightWeight_==1) {
     352      mu_ *= 1000.0;
     353      maxIts2_=23;
     354    } else if (lightWeight_==2) {
     355      maxIts2_=11;
     356    } else {
     357      maxIts2_=23;
     358    }
    141359  }
    142360  //printf("setting mu to %g and doing %d passes\n",mu_,majorIterations_);
     
    144362#ifndef OSI_IDIOT
    145363  double averageInfeas = model_->sumPrimalInfeasibilities()/((double) model_->numberRows());
    146   if (averageInfeas<0.01&&(strategy_&512)!=0)
     364  if ((averageInfeas<0.01&&(strategy_&512)!=0)||(strategy_&8192)!=0)
    147365    crossOver(16+1);
    148366  else
     
    247465  dj=new double[ncols];
    248466  delete [] whenUsed_;
    249   whenUsed=whenUsed_=new int[ncols];
     467  bool oddSlacks=false;
     468  // See if any costed slacks
     469  int numberSlacks=0;
     470  for (i=0;i<ncols;i++) {
     471    if (columnLength[i]==1)
     472      numberSlacks++;
     473  }
     474  if (!numberSlacks) {
     475    whenUsed_=new int[ncols];
     476  } else {
     477    printf("%d slacks\n",numberSlacks);
     478    oddSlacks=true;
     479    int extra = (int) (nrows*sizeof(double)/sizeof(int));
     480    whenUsed_=new int[2*ncols+2*nrows+extra];
     481    int * posSlack = whenUsed_+ncols;
     482    int * negSlack = posSlack+nrows;
     483    int * nextSlack = negSlack + nrows;
     484    for (i=0;i<nrows;i++) {
     485      posSlack[i]=-1;
     486      negSlack[i]=-1;
     487    }
     488    for (i=0;i<ncols;i++)
     489      nextSlack[i]=-1;
     490    for (i=0;i<ncols;i++) {
     491      if (columnLength[i]==1) {
     492        CoinBigIndex j=columnStart[i];
     493        int iRow = row[j];
     494        if (element[j]>0.0) {
     495          if (posSlack[iRow]==-1) {
     496            posSlack[iRow]=i;
     497          } else {
     498            int iCol = posSlack[iRow];
     499            while (nextSlack[iCol]>=0)
     500              iCol = nextSlack[iCol];
     501            nextSlack[iCol]=i;
     502          }
     503        } else {
     504          if (negSlack[iRow]==-1) {
     505            negSlack[iRow]=i;
     506          } else {
     507            int iCol = negSlack[iRow];
     508            while (nextSlack[iCol]>=0)
     509              iCol = nextSlack[iCol];
     510            nextSlack[iCol]=i;
     511          }
     512        }
     513      }
     514    }
     515    // now sort
     516    for (i=0;i<nrows;i++) {
     517      int iCol;
     518      iCol = posSlack[i];
     519      if (iCol>=0) {
     520        CoinBigIndex j=columnStart[iCol];
     521#ifndef NDEBUG
     522        int iRow = row[j];
     523#endif
     524        assert (element[j]>0.0);
     525        assert (iRow==i);
     526        dj[0]= cost[iCol]/element[j];
     527        whenUsed_[0]=iCol;
     528        int n=1;
     529        while (nextSlack[iCol]>=0) {
     530          iCol = nextSlack[iCol];
     531          CoinBigIndex j=columnStart[iCol];
     532#ifndef NDEBUG
     533          int iRow = row[j];
     534#endif
     535          assert (element[j]>0.0);
     536          assert (iRow==i);
     537          dj[n]= cost[iCol]/element[j];
     538          whenUsed_[n++]=iCol;
     539        }
     540        for (j=0;j<n;j++) {
     541          int jCol = whenUsed_[j];
     542          nextSlack[jCol]=-2;
     543        }
     544        CoinSort_2(dj,dj+n,whenUsed_);
     545        // put back
     546        iCol = whenUsed_[0];
     547        posSlack[i]=iCol;
     548        for (j=1;j<n;j++) {
     549          int jCol = whenUsed_[j];
     550          nextSlack[iCol]=jCol;
     551          iCol=jCol;
     552        }
     553      }
     554      iCol = negSlack[i];
     555      if (iCol>=0) {
     556        CoinBigIndex j=columnStart[iCol];
     557#ifndef NDEBUG
     558        int iRow = row[j];
     559#endif
     560        assert (element[j]<0.0);
     561        assert (iRow==i);
     562        dj[0]= -cost[iCol]/element[j];
     563        whenUsed_[0]=iCol;
     564        int n=1;
     565        while (nextSlack[iCol]>=0) {
     566          iCol = nextSlack[iCol];
     567          CoinBigIndex j=columnStart[iCol];
     568#ifndef NDEBUG
     569          int iRow = row[j];
     570#endif
     571          assert (element[j]<0.0);
     572          assert (iRow==i);
     573          dj[n]= -cost[iCol]/element[j];
     574          whenUsed_[n++]=iCol;
     575        }
     576        for (j=0;j<n;j++) {
     577          int jCol = whenUsed_[j];
     578          nextSlack[jCol]=-2;
     579        }
     580        CoinSort_2(dj,dj+n,whenUsed_);
     581        // put back
     582        iCol = whenUsed_[0];
     583        negSlack[i]=iCol;
     584        for (j=1;j<n;j++) {
     585          int jCol = whenUsed_[j];
     586          nextSlack[iCol]=jCol;
     587          iCol=jCol;
     588        }
     589      }
     590    }
     591  }
     592  whenUsed=whenUsed_;
    250593  if (model_->getObjSense()==-1.0) {
    251594    maxmin=-1.0;
     
    385728                       0,mu,drop,
    386729                       maxmin,offset,strategy,djTol,djExit,djFlag);
    387   n=0;
    388   for (i=ordStart;i<ordEnd;i++) {
    389     if (colsol[i]>lower[i]+fixTolerance) {
    390       if (colsol[i]<upper[i]-fixTolerance) {
    391         n++;
    392       } else {
    393         colsol[i]=upper[i];
    394       }
    395       whenUsed[i]=iteration;
    396     } else {
    397       colsol[i]=lower[i];
    398     }
     730  // update whenUsed_
     731  n = cleanIteration(iteration, ordStart,ordEnd,
     732                     colsol,  lower,  upper,
     733                     rowlower, rowupper,
     734                     cost, elemXX, fixTolerance,lastResult.objval,lastResult.infeas);
     735  if ((strategy_&16384)!=0) {
     736    int * posSlack = whenUsed_+ncols;
     737    int * negSlack = posSlack+nrows;
     738    int * nextSlack = negSlack + nrows;
     739    double * rowsol2 = (double *) (nextSlack+ncols);
     740    for (i=0;i<nrows;i++)
     741      rowsol[i] += rowsol2[i];
    399742  }
    400743  if ((logLevel_&1)!=0) {
     
    432775      bestFeasible=lastResult.objval;
    433776    }
     777    if (lastResult.infeas+mu*lastResult.objval<bestWeighted) {
     778      bestWeighted=lastResult.objval+mu*lastResult.objval;
     779    }
    434780    if ((saveStrategy&4096)) strategy |=256;
    435781    if ((saveStrategy&4)!=0&&iteration>2) {
     
    441787      djTol=0.01*djExit;
    442788    }
     789    memcpy(saveSol,colsol,ncols*sizeof(double));
    443790    result=IdiSolve(nrows,ncols,rowsol ,colsol,pi,dj,
    444791                     cost,rowlower,rowupper,
     
    446793                     maxIts,mu,drop,
    447794                     maxmin,offset,strategy,djTol,djExit,djFlag);
    448     n=0;
    449     for (i=ordStart;i<ordEnd;i++) {
    450       if (colsol[i]>lower[i]+fixTolerance) {
    451         if (colsol[i]<upper[i]-fixTolerance) {
    452           n++;
    453         } else {
    454           colsol[i]=upper[i];
    455         }
    456         whenUsed[i]=iteration;
    457       } else {
    458         colsol[i]=lower[i];
    459       }
     795    n = cleanIteration(iteration, ordStart,ordEnd,
     796                       colsol,  lower,  upper,
     797                       rowlower, rowupper,
     798                       cost, elemXX, fixTolerance,result.objval,result.infeas);
     799    if ((strategy_&16384)!=0) {
     800      int * posSlack = whenUsed_+ncols;
     801      int * negSlack = posSlack+nrows;
     802      int * nextSlack = negSlack + nrows;
     803      double * rowsol2 = (double *) (nextSlack+ncols);
     804      for (i=0;i<nrows;i++)
     805        rowsol[i] += rowsol2[i];
    460806    }
    461807    if ((logLevel_&1)!=0) {
     
    473819#endif
    474820    }
    475     if (iteration>50&&n==numberAway&&result.infeas<1.0e-4)
     821    if (iteration>50&&n==numberAway&&result.infeas<1.0e-4) {
     822      printf("infeas small %g\n",result.infeas);
    476823      break; // not much happening
     824    }
    477825    if (lightWeight_==1&&iteration>10&&result.infeas>1.0&&maxIts!=7) {
    478826      if (lastInfeas!=bestInfeas&&CoinMin(result.infeas,lastInfeas)>0.95*bestInfeas)
     
    487835      if ((strategy_&1024)!=0&&mu<1.0e-10)
    488836        result.infeas=firstInfeas*0.8;
    489       if (majorIterations_>=50)
     837      if (majorIterations_>=50||dropEnoughFeasibility_<=0.0)
    490838        result.infeas *= 0.8;
    491839      if (result.infeas>firstInfeas*0.9
     
    506854        memcpy(colsol,saveSol,ncols*sizeof(double));
    507855      } else {
    508         // Save best solution
    509         memcpy(saveSol,colsol,ncols*sizeof(double));
    510856        maxIts=maxIts2;
    511857        checkIteration=0;
    512858        if ((strategy_&1024)!=0) mu *= 1.0e-1;
    513859      }
    514     } else if (result.infeas<bestInfeas) {
    515       // Save best solution
    516       memcpy(saveSol,colsol,ncols*sizeof(double));
     860    } else {
    517861    }
    518862    bestInfeas=CoinMin(bestInfeas,result.infeas);
     
    533877            (result.infeas>lastResult.infeas*0.9
    534878             &&result.weighted>lastResult.weighted
    535              -dropEnoughWeighted_*fabs(lastResult.weighted))) {
     879             -dropEnoughWeighted_*CoinMax(fabs(lastResult.weighted),fabs(result.weighted)))) {
    536880          mu*=changeMu;
    537881          if ((saveStrategy&32)!=0&&result.infeas<reasonableInfeas&&0) {
     
    539883            printf("reasonable infeas now %g\n",reasonableInfeas);
    540884          }
     885          result.weighted=1.0e60;
    541886          nTry=0;
    542887          bestFeasible=1.0e31;
     
    630975    }
    631976    if (iteration>=majorIterations_) {
    632       // If small and not feasible and crash then dive dive dive
    633       if (0&&result.infeas>1.0&&majorIterations_<30&&(maxIts2_==11||maxIts2_==23)) {
    634         maxIts=7;
    635         majorIterations_=100;
     977      // If not feasible and crash then dive dive dive
     978      if (mu>1.0e-12&&result.infeas>1.0&&majorIterations_<40) {
     979        mu=1.0e-30;
     980        majorIterations_=iteration+1;
     981        stopMu=0.0;
    636982      } else {
    637983        if (logLevel>2)
     
    642988  }
    643989  majorIterations_ = saveMajorIterations;
    644   // put back best solution
    645   memcpy(colsol,saveSol,ncols*sizeof(double));
    646990#ifndef OSI_IDIOT
    647991  if (scaled) {
     
    652996    for (icol=0;icol<ncols;icol++) {
    653997      colsol[icol] *= columnScale[icol];
     998      saveSol[icol] *= columnScale[icol];
    654999      dj[icol] /= columnScale[icol];
    6551000    }
     
    7021047  }
    7031048  muAtExit_=mu;
    704   n=0;
    705   for (i=ordStart;i<ordEnd;i++) {
    706     if (colsol[i]>lower[i]+fixTolerance) {
    707       n++;
    708       whenUsed[i]=iteration;
    709     } else {
    710       colsol[i]=lower[i];
    711     }
    712   }
     1049  // For last iteration make as feasible as possible
     1050  if (oddSlacks)
     1051    strategy_ |= 16384;
     1052  // not scaled
     1053  n = cleanIteration(iteration, ordStart,ordEnd,
     1054                     colsol,  lower,  upper,
     1055                     model_->rowLower(), model_->rowUpper(),
     1056                     cost, element, fixTolerance,lastResult.objval,lastResult.infeas);
    7131057#if 0
    714   if ((logLevel&1)==0) {
     1058  if ((logLevel&1)==0||(strategy_&16384)!=0) {
    7151059    printf(
    7161060            "%d - mu %g, infeasibility %g, objective %g, %d interior\n",
     
    7211065  model_->setSumPrimalInfeasibilities(lastResult.infeas);
    7221066#endif
    723   {
     1067  // Put back more feasible solution
     1068  double saveInfeas[]={0.0,0.0};
     1069  for (int iSol=0;iSol<3;iSol++) {
     1070    const double * solution = iSol ? colsol : saveSol;
     1071    if (iSol==2&&saveInfeas[0]<saveInfeas[1]) {
     1072      // put back best solution
     1073      memcpy(colsol,saveSol,ncols*sizeof(double));
     1074    }
    7241075    double large=0.0;
    7251076    int i;
     
    7271078    for (i=0;i<ncols;i++) {
    7281079      CoinBigIndex j;
    729       double value=colsol[i];
     1080      double value=solution[i];
    7301081      for (j=columnStart[i];j<columnStart[i]+columnLength[i];j++) {
    7311082        int irow=row[j];
     
    7441095      }
    7451096    }
     1097    if (iSol<2)
     1098      saveInfeas[iSol]=large;
    7461099    if (logLevel>2)
    7471100      printf("largest infeasibility is %g\n",large);
     
    8411194  double * saveUpper = NULL;
    8421195  double * saveLower = NULL;
     1196  double * saveRowUpper = NULL;
     1197  double * saveRowLower = NULL;
     1198  bool allowInfeasible = (strategy_&8192)!=0;
    8431199  if (addAll<3) {
    8441200    saveUpper = new double [ncols];
     
    8461202    memcpy(saveUpper,upper,ncols*sizeof(double));
    8471203    memcpy(saveLower,lower,ncols*sizeof(double));
     1204    if (allowInfeasible) {
     1205      saveRowUpper = new double [nrows];
     1206      saveRowLower = new double [nrows];
     1207      memcpy(saveRowUpper,rowupper,nrows*sizeof(double));
     1208      memcpy(saveRowLower,rowlower,nrows*sizeof(double));
     1209      double averageInfeas = model_->sumPrimalInfeasibilities()/((double) model_->numberRows());
     1210      fixTolerance = CoinMax(fixTolerance,1.0e-5*averageInfeas);
     1211    }
    8481212  }
    8491213  if (slackStart>=0) {
     
    9291293    }
    9301294    model_->setColumnStatus(i,ClpSimplex::superBasic);
     1295  }
     1296  if ((strategy_&16384)!=0) {
     1297    // put in basis
     1298    int * posSlack = whenUsed_+ncols;
     1299    int * negSlack = posSlack+nrows;
     1300    int * nextSlack = negSlack + nrows;
     1301    for (i=0;i<nrows;i++) {
     1302      int n=0;
     1303      int iCol;
     1304      iCol =posSlack[i];
     1305      if (iCol>=0) {
     1306        if (colsol[iCol]>lower[iCol]+1.0e-8&&
     1307            colsol[iCol]<upper[iCol]+1.0e-8) {
     1308          model_->setColumnStatus(iCol,ClpSimplex::basic);
     1309          n++;
     1310        }
     1311        while (nextSlack[iCol]>=0) {
     1312          iCol = nextSlack[iCol];
     1313          if (colsol[iCol]>lower[iCol]+1.0e-8&&
     1314              colsol[iCol]<upper[iCol]+1.0e-8) {
     1315            model_->setColumnStatus(iCol,ClpSimplex::basic);
     1316            n++;
     1317          }
     1318        }
     1319      }
     1320      iCol =negSlack[i];
     1321      if (iCol>=0) {
     1322        if (colsol[iCol]>lower[iCol]+1.0e-8&&
     1323            colsol[iCol]<upper[iCol]+1.0e-8) {
     1324          model_->setColumnStatus(iCol,ClpSimplex::basic);
     1325          n++;
     1326        }
     1327        while (nextSlack[iCol]>=0) {
     1328          iCol = nextSlack[iCol];
     1329          if (colsol[iCol]>lower[iCol]+1.0e-8&&
     1330              colsol[iCol]<upper[iCol]+1.0e-8) {
     1331            model_->setColumnStatus(iCol,ClpSimplex::basic);
     1332            n++;
     1333          }
     1334        }
     1335      }
     1336      if (n) {
     1337        if (fabs(rowsol[i]-rowlower[i])<fabs(rowsol[i]-rowupper[i]))
     1338          model_->setRowStatus(i,ClpSimplex::atLowerBound);
     1339        else
     1340          model_->setRowStatus(i,ClpSimplex::atUpperBound);
     1341        if (n>1)
     1342          printf("%d basic on row %d!\n",n,i);
     1343      }
     1344    }
    9311345  }
    9321346  double maxmin;
     
    10151429    /*printf("%d in basis\n",ninbas);*/
    10161430  }
     1431  bool wantVector=false;
     1432  if (dynamic_cast< ClpPackedMatrix*>(model_->clpMatrix())) {
     1433    // See if original wanted vector
     1434    ClpPackedMatrix * clpMatrixO = dynamic_cast< ClpPackedMatrix*>(model_->clpMatrix());
     1435    wantVector = clpMatrixO->wantsSpecialColumnCopy();
     1436  }
    10171437  if (addAll<3) {
    10181438    ClpPresolve pinfo;
    10191439    if (presolve) {
     1440      if (allowInfeasible) {
     1441        // fix up so will be feasible
     1442        double * rhs = new double[nrows];
     1443        memset(rhs,0,nrows*sizeof(double));
     1444        model_->clpMatrix()->times(1.0,colsol,rhs);
     1445        double * rowupper = model_->rowUpper();
     1446        double * rowlower= model_->rowLower();
     1447        double sum = 0.0;
     1448        for (i=0;i<nrows;i++) {
     1449          if (rhs[i]>rowupper[i]) {
     1450            sum += rhs[i]-rowupper[i];
     1451            rowupper[i]=rhs[i];
     1452          }
     1453          if (rhs[i]<rowlower[i]) {
     1454            sum += rowlower[i]-rhs[i];
     1455            rowlower[i]=rhs[i];
     1456          }
     1457        }
     1458        printf("sum of infeasibilities %g\n",sum);
     1459        delete [] rhs;
     1460      }
    10201461      saveModel = model_;
     1462      pinfo.setPresolveActions(pinfo.presolveActions()|16384);
    10211463      model_ = pinfo.presolvedModel(*model_,1.0e-8,false,5);
    10221464    }
    10231465    if (model_) {
    1024       model_->primal(1);
     1466      if (!wantVector) {
     1467        model_->primal(1);
     1468      } else {
     1469        ClpMatrixBase * matrix = model_->clpMatrix();
     1470        ClpPackedMatrix * clpMatrix = dynamic_cast< ClpPackedMatrix*>(matrix);
     1471        assert (clpMatrix);
     1472        clpMatrix->makeSpecialColumnCopy();
     1473        model_->primal(1);
     1474        clpMatrix->releaseSpecialColumnCopy();
     1475      }
    10251476      if (presolve) {
    10261477        pinfo.postsolve(true);
     
    10351486      model_ = saveModel;
    10361487      saveModel=NULL;
     1488    }
     1489    if (allowInfeasible) {
     1490      memcpy(model_->rowUpper(),saveRowUpper,nrows*sizeof(double));
     1491      memcpy(model_->rowLower(),saveRowLower,nrows*sizeof(double));
     1492      delete [] saveRowUpper;
     1493      delete [] saveRowLower;
     1494      saveRowUpper = NULL;
     1495      saveRowLower = NULL;
    10371496    }
    10381497    if (addAll<2) {
     
    10581517          }
    10591518        }
     1519        delete [] saveUpper;
     1520        delete [] saveLower;
     1521        saveUpper=NULL;
     1522        saveLower=NULL;
    10601523      }
    10611524      printf("Time so far %g, %d now added from previous iterations\n",
     
    10691532        presolve=0;
    10701533      }
    1071       model_->primal(1);
     1534      if (!wantVector) {
     1535        model_->primal(1);
     1536      } else {
     1537        ClpMatrixBase * matrix = model_->clpMatrix();
     1538        ClpPackedMatrix * clpMatrix = dynamic_cast< ClpPackedMatrix*>(matrix);
     1539        assert (clpMatrix);
     1540        clpMatrix->makeSpecialColumnCopy();
     1541        model_->primal(1);
     1542        clpMatrix->releaseSpecialColumnCopy();
     1543      }
    10721544      if (presolve) {
    10731545        pinfo.postsolve(true);
     
    10851557          }
    10861558        }
     1559        delete [] saveUpper;
     1560        delete [] saveLower;
     1561        saveUpper=NULL;
     1562        saveLower=NULL;
    10871563        printf("Time so far %g, %d now added from previous iterations\n",
    10881564               CoinCpuTime()-startTime,n);
     
    10941570        presolve=0;
    10951571      }
    1096       model_->primal(1);
     1572      if (!wantVector) {
     1573        model_->primal(1);
     1574      } else {
     1575        ClpMatrixBase * matrix = model_->clpMatrix();
     1576        ClpPackedMatrix * clpMatrix = dynamic_cast< ClpPackedMatrix*>(matrix);
     1577        assert (clpMatrix);
     1578        clpMatrix->makeSpecialColumnCopy();
     1579        model_->primal(1);
     1580        clpMatrix->releaseSpecialColumnCopy();
     1581      }
    10971582      if (presolve) {
    10981583        pinfo.postsolve(true);
     
    11431628  baseIts *= 10;
    11441629  maxIts2_ =200+baseIts+5;
     1630  maxIts2_=100;
    11451631  reasonableInfeas_ =((double) nrows)*0.05;
    11461632  lightWeight_=0;
     
    11821668  baseIts *= 10;
    11831669  maxIts2_ =200+baseIts+5;
     1670  maxIts2_=100;
    11841671  reasonableInfeas_ =((double) nrows)*0.05;
    11851672  lightWeight_=0;
Note: See TracChangeset for help on using the changeset viewer.