Changeset 238 for trunk/ClpPresolve.cpp


Ignore:
Timestamp:
Oct 22, 2003 4:45:04 PM (16 years ago)
Author:
forrest
Message:

Messages and stuff

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/ClpPresolve.cpp

    r225 r238  
    4646  nrows_(0),
    4747  nelems_(0),
    48   numberPasses_(5)
     48  numberPasses_(5),
     49  saveFile_("")
    4950{
    5051}
     
    8182                         bool dropNames)
    8283{
    83   ncols_ = si.getNumCols();
    84   nrows_ = si.getNumRows();
    85   nelems_ = si.getNumElements();
     84  // Check matrix
     85  if (!si.clpMatrix()->allElementsInRange(&si,si.getSmallElementValue(),
     86                                          1.0e20))
     87    return NULL;
     88  else
     89    return gutsOfPresolvedModel(&si,feasibilityTolerance,keepIntegers,numberPasses,dropNames);
     90}
     91/* This version of presolve returns a pointer to a new presolved
     92   model and save original data to file.  Returns non-zero if infeasible
     93*/
     94int
     95ClpPresolve::presolvedModel(ClpSimplex &si,std::string fileName,
     96                            double feasibilityTolerance,
     97                            bool keepIntegers,
     98                            int numberPasses)
     99{
     100  // Check matrix
     101  if (!si.clpMatrix()->allElementsInRange(&si,si.getSmallElementValue(),
     102                                          1.0e20))
     103    return 2;
     104  saveFile_=fileName;
     105  si.saveModel(saveFile_.c_str());
     106  ClpSimplex * model = gutsOfPresolvedModel(&si,feasibilityTolerance,keepIntegers,numberPasses,true);
     107  if (model==&si) {
     108    return 0;
     109  } else {
     110    si.restoreModel(saveFile_.c_str());
     111    return 1;
     112  }
     113}
     114
     115// Return pointer to presolved model
     116ClpSimplex *
     117ClpPresolve::model() const
     118{
     119  return presolvedModel_;
     120}
     121// Return pointer to original model
     122ClpSimplex *
     123ClpPresolve::originalModel() const
     124{
     125  return originalModel_;
     126}
     127void
     128ClpPresolve::postsolve(bool updateStatus)
     129{
     130  // Messages
     131  CoinMessages messages = CoinMessage(presolvedModel_->messages().language());
     132  if (!presolvedModel_->isProvenOptimal()) {
     133    presolvedModel_->messageHandler()->message(COIN_PRESOLVE_NONOPTIMAL,
     134                                             messages)
     135                                               <<CoinMessageEol;
     136  }
     137
     138  // this is the size of the original problem
     139  const int ncols0  = ncols_;
     140  const int nrows0  = nrows_;
     141  const CoinBigIndex nelems0 = nelems_;
     142 
     143  // this is the reduced problem
     144  int ncols = presolvedModel_->getNumCols();
     145  int nrows = presolvedModel_->getNumRows();
     146
     147  double * acts=NULL;
     148  double * sol =NULL;
     149  unsigned char * rowstat=NULL;
     150  unsigned char * colstat = NULL;
     151  if (saveFile_=="") {
     152    // reality check
     153    assert(ncols0==originalModel_->getNumCols());
     154    assert(nrows0==originalModel_->getNumRows());
     155    acts = originalModel_->primalRowSolution();
     156    sol  = originalModel_->primalColumnSolution();
     157    if (updateStatus) {
     158      unsigned char *status = originalModel_->statusArray();
     159      rowstat = status + ncols0;
     160      colstat = status;
     161      memcpy(colstat, presolvedModel_->statusArray(), ncols);
     162      memcpy(rowstat, presolvedModel_->statusArray()+ncols, nrows);
     163    }
     164  } else {
     165    // from file
     166    acts = new double[nrows0];
     167    sol  = new double[ncols0];
     168    if (updateStatus) {
     169      unsigned char *status = new unsigned char [nrows0+ncols0];
     170      rowstat = status + ncols0;
     171      colstat = status;
     172      memcpy(colstat, presolvedModel_->statusArray(), ncols);
     173      memcpy(rowstat, presolvedModel_->statusArray()+ncols, nrows);
     174    }
     175  }
     176
     177
     178  CoinPostsolveMatrix prob(presolvedModel_,
     179                       ncols0,
     180                       nrows0,
     181                       nelems0,
     182                       presolvedModel_->getObjSense(),
     183                       // end prepost
     184                       
     185                       sol, acts,
     186                       colstat, rowstat);
     187   
     188  postsolve(prob);
     189
     190  if (saveFile_!="") {
     191    // From file
     192    assert (originalModel_==presolvedModel_);
     193    originalModel_->restoreModel(saveFile_.c_str());
     194    memcpy(originalModel_->primalRowSolution(),acts,nrows0*sizeof(double));
     195    delete [] acts;
     196    memcpy(originalModel_->primalColumnSolution(),sol,ncols0*sizeof(double));
     197    delete [] sol;
     198    if (updateStatus) {
     199      memcpy(originalModel_->statusArray(),colstat,nrows0+ncols0);
     200      delete [] colstat;
     201    }
     202  }
     203  // put back duals
     204  memcpy(originalModel_->dualRowSolution(),prob.rowduals_,
     205         nrows_*sizeof(double));
     206  double maxmin = originalModel_->getObjSense();
     207  if (maxmin<0.0) {
     208    // swap signs
     209    int i;
     210    double * pi = originalModel_->dualRowSolution();
     211    for (i=0;i<nrows_;i++)
     212      pi[i] = -pi[i];
     213  }
     214  // Now check solution
     215  memcpy(originalModel_->dualColumnSolution(),
     216         originalModel_->objective(),ncols_*sizeof(double));
     217  originalModel_->transposeTimes(-1.0,
     218                                 originalModel_->dualRowSolution(),
     219                                 originalModel_->dualColumnSolution());
     220  memset(originalModel_->primalRowSolution(),0,nrows_*sizeof(double));
     221  originalModel_->times(1.0,originalModel_->primalColumnSolution(),
     222                        originalModel_->primalRowSolution());
     223  originalModel_->checkSolution();
     224  // Messages
     225  presolvedModel_->messageHandler()->message(COIN_PRESOLVE_POSTSOLVE,
     226                                            messages)
     227                                              <<originalModel_->objectiveValue()
     228                                              <<originalModel_->sumDualInfeasibilities()
     229                                              <<originalModel_->numberDualInfeasibilities()
     230                                              <<originalModel_->sumPrimalInfeasibilities()
     231                                              <<originalModel_->numberPrimalInfeasibilities()
     232                                               <<CoinMessageEol;
     233 
     234  //originalModel_->objectiveValue_=objectiveValue_;
     235  originalModel_->setNumberIterations(presolvedModel_->numberIterations());
     236  if (!presolvedModel_->status()) {
     237    if (!originalModel_->numberDualInfeasibilities()&&
     238        !originalModel_->numberPrimalInfeasibilities()) {
     239      originalModel_->setProblemStatus( 0);
     240    } else {
     241      originalModel_->setProblemStatus( -1);
     242      presolvedModel_->messageHandler()->message(COIN_PRESOLVE_NEEDS_CLEANING,
     243                                            messages)
     244                                              <<CoinMessageEol;
     245    }
     246  } else {
     247    originalModel_->setProblemStatus( presolvedModel_->status());
     248  }
     249  if (saveFile_!="")
     250    presolvedModel_=NULL;
     251}
     252
     253// return pointer to original columns
     254const int *
     255ClpPresolve::originalColumns() const
     256{
     257  return originalColumn_;
     258}
     259// return pointer to original rows
     260const int *
     261ClpPresolve::originalRows() const
     262{
     263  return originalRow_;
     264}
     265// Set pointer to original model
     266void
     267ClpPresolve::setOriginalModel(ClpSimplex * model)
     268{
     269  originalModel_=model;
     270}
     271#if 0
     272// A lazy way to restrict which transformations are applied
     273// during debugging.
     274static int ATOI(const char *name)
     275{
     276 return true;
     277#if     DEBUG_PRESOLVE || PRESOLVE_SUMMARY
     278  if (getenv(name)) {
     279    int val = atoi(getenv(name));
     280    printf("%s = %d\n", name, val);
     281    return (val);
     282  } else {
     283    if (strcmp(name,"off"))
     284      return (true);
     285    else
     286      return (false);
     287  }
     288#else
     289  return (true);
     290#endif
     291}
     292#endif
     293//#define DEBUG_PRESOLVE 1
     294#if DEBUG_PRESOLVE
     295void check_sol(CoinPresolveMatrix *prob,double tol)
     296{
     297  double *colels        = prob->colels_;
     298  int *hrow             = prob->hrow_;
     299  int *mcstrt           = prob->mcstrt_;
     300  int *hincol           = prob->hincol_;
     301  int *hinrow           = prob->hinrow_;
     302  int ncols             = prob->ncols_;
     303
     304
     305  double * csol = prob->sol_;
     306  double * acts = prob->acts_;
     307  double * clo = prob->clo_;
     308  double * cup = prob->cup_;
     309  int nrows = prob->nrows_;
     310  double * rlo = prob->rlo_;
     311  double * rup = prob->rup_;
     312
     313  int colx;
     314
     315  double * rsol = new double[nrows];
     316  memset(rsol,0,nrows*sizeof(double));
     317
     318  for (colx = 0; colx < ncols; ++colx) {
     319    if (1) {
     320      CoinBigIndex k = mcstrt[colx];
     321      int nx = hincol[colx];
     322      double solutionValue = csol[colx];
     323      for (int i=0; i<nx; ++i) {
     324        int row = hrow[k];
     325        double coeff = colels[k];
     326        k++;
     327        rsol[row] += solutionValue*coeff;
     328      }
     329      if (csol[colx]<clo[colx]-tol) {
     330        printf("low CSOL:  %d  - %g %g %g\n",
     331                   colx, clo[colx], csol[colx], cup[colx]);
     332      } else if (csol[colx]>cup[colx]+tol) {
     333        printf("high CSOL:  %d  - %g %g %g\n",
     334                   colx, clo[colx], csol[colx], cup[colx]);
     335      }
     336    }
     337  }
     338  int rowx;
     339  for (rowx = 0; rowx < nrows; ++rowx) {
     340    if (hinrow[rowx]) {
     341      if (fabs(rsol[rowx]-acts[rowx])>tol)
     342        printf("inacc RSOL:  %d - %g %g (acts_ %g) %g\n",
     343                   rowx,  rlo[rowx], rsol[rowx], acts[rowx], rup[rowx]);
     344      if (rsol[rowx]<rlo[rowx]-tol) {
     345        printf("low RSOL:  %d - %g %g %g\n",
     346                   rowx,  rlo[rowx], rsol[rowx], rup[rowx]);
     347      } else if (rsol[rowx]>rup[rowx]+tol ) {
     348        printf("high RSOL:  %d - %g %g %g\n",
     349                   rowx,  rlo[rowx], rsol[rowx], rup[rowx]);
     350      }
     351    }
     352  }
     353  delete [] rsol;
     354}
     355#endif
     356// This is the presolve loop.
     357// It is a separate virtual function so that it can be easily
     358// customized by subclassing CoinPresolve.
     359const CoinPresolveAction *ClpPresolve::presolve(CoinPresolveMatrix *prob)
     360{
     361  paction_ = 0;
     362
     363  prob->status_=0; // say feasible
     364
     365  paction_ = make_fixed(prob, paction_);
     366  // if integers then switch off dual stuff
     367  // later just do individually
     368  bool doDualStuff = (presolvedModel_->integerInformation()==NULL);
     369
     370#if     CHECK_CONSISTENCY
     371  presolve_links_ok(prob->rlink_, prob->mrstrt_, prob->hinrow_, prob->nrows_);
     372#endif
     373
     374  if (!prob->status_) {
     375#if 0
     376    const bool slackd = ATOI("SLACKD")!=0;
     377    //const bool forcing = ATOI("FORCING")!=0;
     378    const bool doubleton = ATOI("DOUBLETON")!=0;
     379    const bool forcing = ATOI("off")!=0;
     380    const bool ifree = ATOI("off")!=0;
     381    const bool zerocost = ATOI("off")!=0;
     382    const bool dupcol = ATOI("off")!=0;
     383    const bool duprow = ATOI("off")!=0;
     384    const bool dual = ATOI("off")!=0;
     385#else
     386    // normal
     387#if 0
     388    const bool slackd = true;
     389    const bool doubleton = false;
     390    const bool tripleton = true;
     391    const bool forcing = false;
     392    const bool ifree = true;
     393    const bool zerocost = true;
     394    const bool dupcol = false;
     395    const bool duprow = false;
     396    const bool dual = doDualStuff;
     397#else
     398    const bool slackd = true;
     399    const bool doubleton = true;
     400    const bool tripleton = true;
     401    const bool forcing = true;
     402    const bool ifree = true;
     403    const bool zerocost = true;
     404    const bool dupcol = true;
     405    const bool duprow = true;
     406    const bool dual = doDualStuff;
     407#endif
     408#endif
     409   
     410    // some things are expensive so just do once (normally)
     411
     412    int i;
     413    // say look at all
     414    if (!prob->anyProhibited()) {
     415      for (i=0;i<nrows_;i++)
     416        prob->rowsToDo_[i]=i;
     417      prob->numberRowsToDo_=nrows_;
     418      for (i=0;i<ncols_;i++)
     419        prob->colsToDo_[i]=i;
     420      prob->numberColsToDo_=ncols_;
     421    } else {
     422      // some stuff must be left alone
     423      prob->numberRowsToDo_=0;
     424      for (i=0;i<nrows_;i++)
     425        if (!prob->rowProhibited(i))
     426            prob->rowsToDo_[prob->numberRowsToDo_++]=i;
     427      prob->numberColsToDo_=0;
     428      for (i=0;i<ncols_;i++)
     429        if (!prob->colProhibited(i))
     430            prob->colsToDo_[prob->numberColsToDo_++]=i;
     431    }
     432
     433
     434    int iLoop;
     435#if     DEBUG_PRESOLVE
     436    check_sol(prob,1.0e0);
     437#endif
     438
     439    // Check number rows dropped
     440    int lastDropped=0;
     441    for (iLoop=0;iLoop<numberPasses_;iLoop++) {
     442#ifdef PRESOLVE_SUMMARY
     443      printf("Starting major pass %d\n",iLoop+1);
     444#endif
     445      const CoinPresolveAction * const paction0 = paction_;
     446      // look for substitutions with no fill
     447      int fill_level=2;
     448      //fill_level=10;
     449      //printf("** fill_level == 10 !\n");
     450      int whichPass=0;
     451      while (1) {
     452        whichPass++;
     453        const CoinPresolveAction * const paction1 = paction_;
     454
     455        if (slackd) {
     456          bool notFinished = true;
     457          while (notFinished)
     458            paction_ = slack_doubleton_action::presolve(prob, paction_,
     459                                                        notFinished);
     460          if (prob->status_)
     461            break;
     462        }
     463
     464        if (doubleton) {
     465          paction_ = doubleton_action::presolve(prob, paction_);
     466          if (prob->status_)
     467            break;
     468        }
     469
     470        if (tripleton) {
     471          paction_ = tripleton_action::presolve(prob, paction_);
     472          if (prob->status_)
     473            break;
     474        }
     475
     476        if (zerocost) {
     477          paction_ = do_tighten_action::presolve(prob, paction_);
     478          if (prob->status_)
     479            break;
     480        }
     481
     482        if (forcing) {
     483          paction_ = forcing_constraint_action::presolve(prob, paction_);
     484          if (prob->status_)
     485            break;
     486        }
     487
     488        if (ifree) {
     489          paction_ = implied_free_action::presolve(prob, paction_,fill_level);
     490          if (prob->status_)
     491            break;
     492        }
     493
     494#if     DEBUG_PRESOLVE
     495        check_sol(prob,1.0e0);
     496#endif
     497
     498#if     CHECK_CONSISTENCY
     499        presolve_links_ok(prob->rlink_, prob->mrstrt_, prob->hinrow_,
     500                          prob->nrows_);
     501#endif
     502
     503#if     DEBUG_PRESOLVE
     504        presolve_no_zeros(prob->mcstrt_, prob->colels_, prob->hincol_,
     505                          prob->ncols_);
     506#endif
     507#if     CHECK_CONSISTENCY
     508        prob->consistent();
     509#endif
     510
     511         
     512        // set up for next pass
     513        // later do faster if many changes i.e. memset and memcpy
     514        prob->numberRowsToDo_ = prob->numberNextRowsToDo_;
     515        int kcheck;
     516        bool found=false;
     517        kcheck=-1;
     518        for (i=0;i<prob->numberNextRowsToDo_;i++) {
     519          int index = prob->nextRowsToDo_[i];
     520          prob->unsetRowChanged(index);
     521          prob->rowsToDo_[i] = index;
     522          if (index==kcheck) {
     523            printf("row %d on list after pass %d\n",kcheck,
     524                   whichPass);
     525            found=true;
     526          }
     527        }
     528        if (!found&&kcheck>=0)
     529          prob->rowsToDo_[prob->numberRowsToDo_++]=kcheck;
     530        prob->numberNextRowsToDo_=0;
     531        prob->numberColsToDo_ = prob->numberNextColsToDo_;
     532        kcheck=-1;
     533        found=false;
     534        for (i=0;i<prob->numberNextColsToDo_;i++) {
     535          int index = prob->nextColsToDo_[i];
     536          prob->unsetColChanged(index);
     537          prob->colsToDo_[i] = index;
     538          if (index==kcheck) {
     539            printf("col %d on list after pass %d\n",kcheck,
     540                   whichPass);
     541            found=true;
     542          }
     543        }
     544        if (!found&&kcheck>=0)
     545          prob->colsToDo_[prob->numberColsToDo_++]=kcheck;
     546        prob->numberNextColsToDo_=0;
     547        if (paction_ == paction1&&fill_level>0)
     548          break;
     549      }
     550      // say look at all
     551      int i;
     552      if (!prob->anyProhibited()) {
     553        for (i=0;i<nrows_;i++)
     554          prob->rowsToDo_[i]=i;
     555        prob->numberRowsToDo_=nrows_;
     556        for (i=0;i<ncols_;i++)
     557          prob->colsToDo_[i]=i;
     558        prob->numberColsToDo_=ncols_;
     559      } else {
     560        // some stuff must be left alone
     561        prob->numberRowsToDo_=0;
     562        for (i=0;i<nrows_;i++)
     563          if (!prob->rowProhibited(i))
     564            prob->rowsToDo_[prob->numberRowsToDo_++]=i;
     565        prob->numberColsToDo_=0;
     566        for (i=0;i<ncols_;i++)
     567          if (!prob->colProhibited(i))
     568            prob->colsToDo_[prob->numberColsToDo_++]=i;
     569      }
     570      // now expensive things
     571      // this caused world.mps to run into numerical difficulties
     572#ifdef PRESOLVE_SUMMARY
     573      printf("Starting expensive\n");
     574#endif
     575
     576      if (dual) {
     577        int itry;
     578        for (itry=0;itry<5;itry++) {
     579          const CoinPresolveAction * const paction2 = paction_;
     580          paction_ = remove_dual_action::presolve(prob, paction_);
     581          if (prob->status_)
     582            break;
     583          if (ifree) {
     584            int fill_level=0; // switches off substitution
     585            paction_ = implied_free_action::presolve(prob, paction_,fill_level);
     586            if (prob->status_)
     587              break;
     588          }
     589          if (paction_ == paction2)
     590            break;
     591        }
     592      }
     593#if     DEBUG_PRESOLVE
     594      check_sol(prob,1.0e0);
     595#endif
     596      if (dupcol) {
     597        paction_ = dupcol_action::presolve(prob, paction_);
     598        if (prob->status_)
     599          break;
     600      }
     601#if     DEBUG_PRESOLVE
     602        check_sol(prob,1.0e0);
     603#endif
     604     
     605      if (duprow) {
     606        paction_ = duprow_action::presolve(prob, paction_);
     607        if (prob->status_)
     608          break;
     609      }
     610#if     DEBUG_PRESOLVE
     611      check_sol(prob,1.0e0);
     612#endif
     613      {
     614        int * hinrow = prob->hinrow_;
     615        int numberDropped=0;
     616        for (i=0;i<nrows_;i++)
     617          if (!hinrow[i])
     618            numberDropped++;
     619        //printf("%d rows dropped after pass %d\n",numberDropped,
     620        //     iLoop+1);
     621        if (numberDropped==lastDropped)
     622          break;
     623        else
     624          lastDropped = numberDropped;
     625      }
     626      if (paction_ == paction0)
     627        break;
     628         
     629    }
     630  }
     631  if (!prob->status_) {
     632    paction_ = drop_zero_coefficients(prob, paction_);
     633#if     DEBUG_PRESOLVE
     634        check_sol(prob,1.0e0);
     635#endif
     636
     637    paction_ = drop_empty_cols_action::presolve(prob, paction_);
     638    paction_ = drop_empty_rows_action::presolve(prob, paction_);
     639#if     DEBUG_PRESOLVE
     640        check_sol(prob,1.0e0);
     641#endif
     642  }
     643 
     644  // Messages
     645  CoinMessages messages = CoinMessage(prob->messages().language());
     646  if (prob->status_) {
     647    if (prob->status_==1)
     648          prob->messageHandler()->message(COIN_PRESOLVE_INFEAS,
     649                                             messages)
     650                                               <<prob->feasibilityTolerance_
     651                                               <<CoinMessageEol;
     652    else if (prob->status_==2)
     653          prob->messageHandler()->message(COIN_PRESOLVE_UNBOUND,
     654                                             messages)
     655                                               <<CoinMessageEol;
     656    else
     657          prob->messageHandler()->message(COIN_PRESOLVE_INFEASUNBOUND,
     658                                             messages)
     659                                               <<CoinMessageEol;
     660    // get rid of data
     661    gutsOfDestroy();
     662  }
     663  return (paction_);
     664}
     665
     666void check_djs(CoinPostsolveMatrix *prob);
     667
     668
     669// We could have implemented this by having each postsolve routine
     670// directly call the next one, but this may make it easier to add debugging checks.
     671void ClpPresolve::postsolve(CoinPostsolveMatrix &prob)
     672{
     673  {
     674    // Check activities
     675    double *colels      = prob.colels_;
     676    int *hrow           = prob.hrow_;
     677    int *mcstrt         = prob.mcstrt_;
     678    int *hincol         = prob.hincol_;
     679    int *link           = prob.link_;
     680    int ncols           = prob.ncols_;
     681
     682    char *cdone = prob.cdone_;
     683
     684    double * csol = prob.sol_;
     685    int nrows = prob.nrows_;
     686
     687    int colx;
     688
     689    double * rsol = prob.acts_;
     690    memset(rsol,0,nrows*sizeof(double));
     691
     692    for (colx = 0; colx < ncols; ++colx) {
     693      if (cdone[colx]) {
     694        CoinBigIndex k = mcstrt[colx];
     695        int nx = hincol[colx];
     696        double solutionValue = csol[colx];
     697        for (int i=0; i<nx; ++i) {
     698          int row = hrow[k];
     699          double coeff = colels[k];
     700          k = link[k];
     701          rsol[row] += solutionValue*coeff;
     702        }
     703      }
     704    }
     705  }
     706  const CoinPresolveAction *paction = paction_;
     707
     708  if (prob.colstat_)
     709    prob.check_nbasic();
     710 
     711#if     DEBUG_PRESOLVE
     712  check_djs(&prob);
     713#endif
     714 
     715 
     716  while (paction) {
     717#if     DEBUG_PRESOLVE
     718    printf("POSTSOLVING %s\n", paction->name());
     719#endif
     720
     721    paction->postsolve(&prob);
     722   
     723#if     DEBUG_PRESOLVE
     724    if (prob.colstat_)
     725      prob.check_nbasic();
     726#endif
     727    paction = paction->next;
     728#if     DEBUG_PRESOLVE
     729    check_djs(&prob);
     730#endif
     731  }   
     732 
     733#if     0 && DEBUG_PRESOLVE
     734  for (i=0; i<ncols0; i++) {
     735    if (!cdone[i]) {
     736      printf("!cdone[%d]\n", i);
     737      abort();
     738    }
     739  }
     740 
     741  for (i=0; i<nrows0; i++) {
     742    if (!rdone[i]) {
     743      printf("!rdone[%d]\n", i);
     744      abort();
     745    }
     746  }
     747 
     748 
     749  for (i=0; i<ncols0; i++) {
     750    if (sol[i] < -1e10 || sol[i] > 1e10)
     751      printf("!!!%d %g\n", i, sol[i]);
     752   
     753  }
     754 
     755 
     756#endif
     757 
     758#if     0 && DEBUG_PRESOLVE
     759  // debug check:  make sure we ended up with same original matrix
     760  {
     761    int identical = 1;
     762   
     763    for (int i=0; i<ncols0; i++) {
     764      PRESOLVEASSERT(hincol[i] == &prob->mcstrt0[i+1] - &prob->mcstrt0[i]);
     765      CoinBigIndex kcs0 = &prob->mcstrt0[i];
     766      CoinBigIndex kcs = mcstrt[i];
     767      int n = hincol[i];
     768      for (int k=0; k<n; k++) {
     769        CoinBigIndex k1 = presolve_find_row1(&prob->hrow0[kcs0+k], kcs, kcs+n, hrow);
     770
     771        if (k1 == kcs+n) {
     772          printf("ROW %d NOT IN COL %d\n", &prob->hrow0[kcs0+k], i);
     773          abort();
     774        }
     775
     776        if (colels[k1] != &prob->dels0[kcs0+k])
     777          printf("BAD COLEL[%d %d %d]:  %g\n",
     778                 k1, i, &prob->hrow0[kcs0+k], colels[k1] - &prob->dels0[kcs0+k]);
     779
     780        if (kcs0+k != k1)
     781          identical=0;
     782      }
     783    }
     784    printf("identical? %d\n", identical);
     785  }
     786#endif
     787}
     788
     789
     790
     791
     792
     793
     794
     795
     796#if     DEBUG_PRESOLVE
     797void check_djs(CoinPostsolveMatrix *prob)
     798{
     799  //return;
     800  double *colels        = prob->colels_;
     801  int *hrow             = prob->hrow_;
     802  int *mcstrt           = prob->mcstrt_;
     803  int *hincol           = prob->hincol_;
     804  int *link             = prob->link_;
     805  int ncols             = prob->ncols_;
     806
     807  double *dcost = prob->cost_;
     808
     809  double *rcosts        = prob->rcosts_;
     810
     811  double *rowduals = prob->rowduals_;
     812
     813  const double maxmin   = prob->maxmin_;
     814
     815  char *cdone   = prob->cdone_;
     816
     817  double * csol = prob->sol_;
     818  double * clo = prob->clo_;
     819  double * cup = prob->cup_;
     820  int nrows = prob->nrows_;
     821  double * rlo = prob->rlo_;
     822  double * rup = prob->rup_;
     823  char *rdone   = prob->rdone_;
     824
     825  int colx;
     826
     827  double * rsol = new double[nrows];
     828  memset(rsol,0,nrows*sizeof(double));
     829
     830  for (colx = 0; colx < ncols; ++colx) {
     831    if (cdone[colx]) {
     832      CoinBigIndex k = mcstrt[colx];
     833      int nx = hincol[colx];
     834      double dj = maxmin * dcost[colx];
     835      double solutionValue = csol[colx];
     836      for (int i=0; i<nx; ++i) {
     837        int row = hrow[k];
     838        double coeff = colels[k];
     839        k = link[k];
     840        dj -= rowduals[row] * coeff;
     841        rsol[row] += solutionValue*coeff;
     842      }
     843      if (! (fabs(rcosts[colx] - dj) < 1.0e-4))
     844        printf("BAD DJ:  %d %g %g\n",
     845               colx, rcosts[colx], dj);
     846      if (cup[colx]-clo[colx]>1.0e-6) {
     847        if (csol[colx]<clo[colx]+1.0e-6) {
     848          if (dj <-1.0e-6)
     849            printf("neg DJ:  %d %g  - %g %g %g\n",
     850                   colx, dj, clo[colx], csol[colx], cup[colx]);
     851        } else if (csol[colx]>cup[colx]-1.0e-6) {
     852          if (dj > 1.0e-6)
     853            printf("pos DJ:  %d %g  - %g %g %g\n",
     854                   colx, dj, clo[colx], csol[colx], cup[colx]);
     855        } else {
     856          if (fabs(dj) >1.0e-6)
     857            printf("nonzero DJ:  %d %g  - %g %g %g\n",
     858                   colx, dj, clo[colx], csol[colx], cup[colx]);
     859        }
     860      }
     861    }
     862  }
     863  int rowx;
     864  for (rowx = 0; rowx < nrows; ++rowx) {
     865    if (rdone[rowx]) {
     866      if (rup[rowx]-rlo[rowx]>1.0e-6) {
     867        double dj = rowduals[rowx];
     868        if (rsol[rowx]<rlo[rowx]+1.0e-6) {
     869          if (dj <-1.0e-5)
     870            printf("neg rDJ:  %d %g  - %g %g %g\n",
     871                   rowx, dj, rlo[rowx], rsol[rowx], rup[rowx]);
     872        } else if (rsol[rowx]>rup[rowx]-1.0e-6) {
     873          if (dj > 1.0e-5)
     874            printf("pos rDJ:  %d %g  - %g %g %g\n",
     875                   rowx, dj, rlo[rowx], rsol[rowx], rup[rowx]);
     876        } else {
     877          if (fabs(dj) >1.0e-5)
     878            printf("nonzero rDJ:  %d %g  - %g %g %g\n",
     879                   rowx, dj, rlo[rowx], rsol[rowx], rup[rowx]);
     880        }
     881      }
     882    }
     883  }
     884  delete [] rsol;
     885}
     886#endif
     887
     888static inline double getTolerance(const ClpSimplex  *si, ClpDblParam key)
     889{
     890  double tol;
     891  if (! si->getDblParam(key, tol)) {
     892    CoinPresolveAction::throwCoinError("getDblParam failed",
     893                                      "CoinPrePostsolveMatrix::CoinPrePostsolveMatrix");
     894  }
     895  return (tol);
     896}
     897
     898
     899// Assumptions:
     900// 1. nrows>=m.getNumRows()
     901// 2. ncols>=m.getNumCols()
     902//
     903// In presolve, these values are equal.
     904// In postsolve, they may be inequal, since the reduced problem
     905// may be smaller, but we need room for the large problem.
     906// ncols may be larger than si.getNumCols() in postsolve,
     907// this at that point si will be the reduced problem,
     908// but we need to reserve enough space for the original problem.
     909CoinPrePostsolveMatrix::CoinPrePostsolveMatrix(const ClpSimplex * si,
     910                                             int ncols_in,
     911                                             int nrows_in,
     912                                             CoinBigIndex nelems_in) :
     913  ncols_(si->getNumCols()),
     914  ncols0_(ncols_in),
     915  nelems_(si->getNumElements()),
     916
     917  mcstrt_(new CoinBigIndex[ncols_in+1]),
     918  hincol_(new int[ncols_in+1]),
     919  hrow_  (new int   [2*nelems_in]),
     920  colels_(new double[2*nelems_in]),
     921
     922  cost_(new double[ncols_in]),
     923  clo_(new double[ncols_in]),
     924  cup_(new double[ncols_in]),
     925  rlo_(new double[nrows_in]),
     926  rup_(new double[nrows_in]),
     927  originalColumn_(new int[ncols_in]),
     928  originalRow_(new int[nrows_in]),
     929
     930  ztolzb_(getTolerance(si, ClpPrimalTolerance)),
     931  ztoldj_(getTolerance(si, ClpDualTolerance)),
     932
     933  maxmin_(si->getObjSense())
     934
     935{
     936  si->getDblParam(ClpObjOffset,originalOffset_);
     937  int ncols = si->getNumCols();
     938  int nrows = si->getNumRows();
     939
     940  ClpDisjointCopyN(si->getColLower(), ncols, clo_);
     941  ClpDisjointCopyN(si->getColUpper(), ncols, cup_);
     942  ClpDisjointCopyN(si->getObjCoefficients(), ncols, cost_);
     943  ClpDisjointCopyN(si->getRowLower(), nrows,  rlo_);
     944  ClpDisjointCopyN(si->getRowUpper(), nrows,  rup_);
     945  int i;
     946  for (i=0;i<ncols_in;i++)
     947    originalColumn_[i]=i;
     948  for (i=0;i<nrows_in;i++)
     949    originalRow_[i]=i;
     950  sol_=NULL;
     951  rowduals_=NULL;
     952  acts_=NULL;
     953
     954  rcosts_=NULL;
     955  colstat_=NULL;
     956  rowstat_=NULL;
     957}
     958
     959// I am not familiar enough with CoinPackedMatrix to be confident
     960// that I will implement a row-ordered version of toColumnOrderedGapFree
     961// properly.
     962static bool isGapFree(const CoinPackedMatrix& matrix)
     963{
     964  const CoinBigIndex * start = matrix.getVectorStarts();
     965  const int * length = matrix.getVectorLengths();
     966  int i;
     967  for (i = matrix.getSizeVectorLengths() - 1; i >= 0; --i) {
     968    if (start[i+1] - start[i] != length[i])
     969      break;
     970  }
     971  return (! (i >= 0));
     972}
     973#if     DEBUG_PRESOLVE
     974static void matrix_bounds_ok(const double *lo, const double *up, int n)
     975{
     976  int i;
     977  for (i=0; i<n; i++) {
     978    PRESOLVEASSERT(lo[i] <= up[i]);
     979    PRESOLVEASSERT(lo[i] < PRESOLVE_INF);
     980    PRESOLVEASSERT(-PRESOLVE_INF < up[i]);
     981  }
     982}
     983#endif
     984CoinPresolveMatrix::CoinPresolveMatrix(int ncols0_in,
     985                                     double maxmin_,
     986                                     // end prepost members
     987
     988                                     ClpSimplex * si,
     989
     990                                     // rowrep
     991                                     int nrows_in,
     992                                     CoinBigIndex nelems_in,
     993                               bool doStatus,
     994                               double nonLinearValue) :
     995
     996  CoinPrePostsolveMatrix(si,
     997                        ncols0_in, nrows_in, nelems_in),
     998  clink_(new presolvehlink[ncols0_in+1]),
     999  rlink_(new presolvehlink[nrows_in+1]),
     1000
     1001  dobias_(0.0),
     1002
     1003  nrows_(si->getNumRows()),
     1004
     1005  // temporary init
     1006  mrstrt_(new CoinBigIndex[nrows_in+1]),
     1007  hinrow_(new int[nrows_in+1]),
     1008  rowels_(new double[2*nelems_in]),
     1009  hcol_(new int[2*nelems_in]),
     1010  integerType_(new char[ncols0_in]),
     1011  feasibilityTolerance_(0.0),
     1012  status_(-1),
     1013  rowsToDo_(new int [nrows_in]),
     1014  numberRowsToDo_(0),
     1015  nextRowsToDo_(new int[nrows_in]),
     1016  numberNextRowsToDo_(0),
     1017  colsToDo_(new int [ncols0_in]),
     1018  numberColsToDo_(0),
     1019  nextColsToDo_(new int[ncols0_in]),
     1020  numberNextColsToDo_(0)
     1021
     1022{
     1023  const int bufsize = 2*nelems_in;
     1024
     1025  // Set up change bits etc
     1026  rowChanged_ = new unsigned char[nrows_];
     1027  memset(rowChanged_,0,nrows_);
     1028  colChanged_ = new unsigned char[ncols_];
     1029  memset(colChanged_,0,ncols_);
     1030  CoinPackedMatrix * m = si->matrix();
     1031
     1032  // The coefficient matrix is a big hunk of stuff.
     1033  // Do the copy here to try to avoid running out of memory.
     1034
     1035  const CoinBigIndex * start = m->getVectorStarts();
     1036  const int * length = m->getVectorLengths();
     1037  const int * row = m->getIndices();
     1038  const double * element = m->getElements();
     1039  int icol,nel=0;
     1040  mcstrt_[0]=0;
     1041  for (icol=0;icol<ncols_;icol++) {
     1042    int j;
     1043    for (j=start[icol];j<start[icol]+length[icol];j++) {
     1044      hrow_[nel]=row[j];
     1045      colels_[nel++]=element[j];
     1046    }
     1047    mcstrt_[icol+1]=nel;
     1048  }
     1049  assert(mcstrt_[ncols_] == nelems_);
     1050  ClpDisjointCopyN(m->getVectorLengths(),ncols_,  hincol_);
     1051
     1052  // same thing for row rep
     1053  m = new CoinPackedMatrix();
     1054  m->reverseOrderedCopyOf(*si->matrix());
     1055  m->removeGaps();
     1056
     1057
     1058  ClpDisjointCopyN(m->getVectorStarts(),  nrows_,  mrstrt_);
     1059  mrstrt_[nrows_] = nelems_;
     1060  ClpDisjointCopyN(m->getVectorLengths(), nrows_,  hinrow_);
     1061  ClpDisjointCopyN(m->getIndices(),       nelems_, hcol_);
     1062  ClpDisjointCopyN(m->getElements(),      nelems_, rowels_);
     1063
     1064  delete m;
     1065  if (si->integerInformation()) {
     1066    memcpy(integerType_,si->integerInformation(),ncols_*sizeof(char));
     1067  } else {
     1068    ClpFillN<char>(integerType_, ncols_, 0);
     1069  }
     1070
     1071  // Set up prohibited bits if needed
     1072  if (nonLinearValue) {
     1073    anyProhibited_ = true;
     1074    for (icol=0;icol<ncols_;icol++) {
     1075      int j;
     1076      bool nonLinearColumn = false;
     1077      if (cost_[icol]==nonLinearValue)
     1078        nonLinearColumn=true;
     1079      for (j=mcstrt_[icol];j<mcstrt_[icol+1];j++) {
     1080        if (colels_[j]==nonLinearValue) {
     1081          nonLinearColumn=true;
     1082          setRowProhibited(hrow_[j]);
     1083        }
     1084      }
     1085      if (nonLinearColumn)
     1086        setColProhibited(icol);
     1087    }
     1088  } else {
     1089    anyProhibited_ = false;
     1090  }
     1091
     1092  if (doStatus) {
     1093    // allow for status and solution
     1094    sol_ = new double[ncols_];
     1095    memcpy(sol_,si->primalColumnSolution(),ncols_*sizeof(double));;
     1096    acts_ = new double [nrows_];
     1097    memcpy(acts_,si->primalRowSolution(),nrows_*sizeof(double));
     1098    if (!si->statusArray())
     1099      si->createStatus();
     1100    colstat_ = new unsigned char [nrows_+ncols_];
     1101    memcpy(colstat_,si->statusArray(),
     1102           (nrows_+ncols_)*sizeof(unsigned char));
     1103    rowstat_ = colstat_+ncols_;
     1104  }
     1105
     1106  // the original model's fields are now unneeded - free them
     1107 
     1108  si->resize(0,0);
     1109
     1110#if     DEBUG_PRESOLVE
     1111  matrix_bounds_ok(rlo_, rup_, nrows_);
     1112  matrix_bounds_ok(clo_, cup_, ncols_);
     1113#endif
     1114
     1115#if 0
     1116  for (i=0; i<nrows; ++i)
     1117    printf("NR: %6d\n", hinrow[i]);
     1118  for (int i=0; i<ncols; ++i)
     1119    printf("NC: %6d\n", hincol[i]);
     1120#endif
     1121
     1122  presolve_make_memlists(mcstrt_, hincol_, clink_, ncols_);
     1123  presolve_make_memlists(mrstrt_, hinrow_, rlink_, nrows_);
     1124
     1125  // this allows last col/row to expand up to bufsize-1 (22);
     1126  // this must come after the calls to presolve_prefix
     1127  mcstrt_[ncols_] = bufsize-1;
     1128  mrstrt_[nrows_] = bufsize-1;
     1129
     1130#if     CHECK_CONSISTENCY
     1131  consistent(false);
     1132#endif
     1133}
     1134
     1135void CoinPresolveMatrix::update_model(ClpSimplex * si,
     1136                                     int nrows0,
     1137                                     int ncols0,
     1138                                     CoinBigIndex nelems0)
     1139{
     1140  si->loadProblem(ncols_, nrows_, mcstrt_, hrow_, colels_, hincol_,
     1141                 clo_, cup_, cost_, rlo_, rup_);
     1142
     1143  delete [] si->integerInformation();
     1144  int numberIntegers=0;
     1145  for (int i=0; i<ncols_; i++) {
     1146    if (integerType_[i])
     1147      numberIntegers++;
     1148  }
     1149  if (numberIntegers)
     1150    si->copyInIntegerInformation(integerType_);
     1151  else
     1152    si->copyInIntegerInformation(NULL);
     1153
     1154#if     PRESOLVE_SUMMARY
     1155  printf("NEW NCOL/NROW/NELS:  %d(-%d) %d(-%d) %d(-%d)\n",
     1156         ncols_, ncols0-ncols_,
     1157         nrows_, nrows0-nrows_,
     1158         si->getNumElements(), nelems0-si->getNumElements());
     1159#endif
     1160  si->setDblParam(ClpObjOffset,originalOffset_-dobias_);
     1161
     1162}
     1163
     1164
     1165
     1166
     1167
     1168
     1169
     1170
     1171
     1172
     1173
     1174////////////////  POSTSOLVE
     1175
     1176CoinPostsolveMatrix::CoinPostsolveMatrix(ClpSimplex*  si,
     1177                                       int ncols0_in,
     1178                                       int nrows0_in,
     1179                                       CoinBigIndex nelems0,
     1180                                   
     1181                                       double maxmin_,
     1182                                       // end prepost members
     1183
     1184                                       double *sol_in,
     1185                                       double *acts_in,
     1186
     1187                                       unsigned char *colstat_in,
     1188                                       unsigned char *rowstat_in) :
     1189  CoinPrePostsolveMatrix(si,
     1190                        ncols0_in, nrows0_in, nelems0),
     1191
     1192  free_list_(0),
     1193  // link, free_list, maxlink
     1194  maxlink_(2*nelems0),
     1195  link_(new int[/*maxlink*/ 2*nelems0]),
     1196     
     1197  cdone_(new char[ncols0_]),
     1198  rdone_(new char[nrows0_in]),
     1199
     1200  nrows_(si->getNumRows()),
     1201  nrows0_(nrows0_in)
     1202{
     1203
     1204  sol_=sol_in;
     1205  rowduals_=NULL;
     1206  acts_=acts_in;
     1207
     1208  rcosts_=NULL;
     1209  colstat_=colstat_in;
     1210  rowstat_=rowstat_in;
     1211
     1212  // this is the *reduced* model, which is probably smaller
     1213  int ncols1 = si->getNumCols();
     1214  int nrows1 = si->getNumRows();
     1215
     1216  const CoinPackedMatrix * m = si->matrix();
     1217
     1218  if (! isGapFree(*m)) {
     1219    CoinPresolveAction::throwCoinError("Matrix not gap free",
     1220                                      "CoinPostsolveMatrix");
     1221  }
     1222
     1223  const CoinBigIndex nelemsr = m->getNumElements();
     1224
     1225  ClpDisjointCopyN(m->getVectorStarts(), ncols1, mcstrt_);
     1226  mcstrt_[ncols_] = nelems0;    // ??
     1227  ClpDisjointCopyN(m->getVectorLengths(),ncols1,  hincol_);
     1228  ClpDisjointCopyN(m->getIndices(),      nelemsr, hrow_);
     1229  ClpDisjointCopyN(m->getElements(),     nelemsr, colels_);
     1230
     1231
     1232#if     0 && DEBUG_PRESOLVE
     1233  presolve_check_costs(model, &colcopy);
     1234#endif
     1235
     1236  // This determines the size of the data structure that contains
     1237  // the matrix being postsolved.  Links are taken from the free_list
     1238  // to recreate matrix entries that were presolved away,
     1239  // and links are added to the free_list when entries created during
     1240  // presolve are discarded.  There is never a need to gc this list.
     1241  // Naturally, it should contain
     1242  // exactly nelems0 entries "in use" when postsolving is done,
     1243  // but I don't know whether the matrix could temporarily get
     1244  // larger during postsolving.  Substitution into more than two
     1245  // rows could do that, in principle.  I am being very conservative
     1246  // here by reserving much more than the amount of space I probably need.
     1247  // If this guess is wrong, check_free_list may be called.
     1248  //  int bufsize = 2*nelems0;
     1249
     1250  memset(cdone_, -1, ncols0_);
     1251  memset(rdone_, -1, nrows0_);
     1252
     1253  rowduals_ = new double[nrows0_];
     1254  ClpDisjointCopyN(si->getRowPrice(), nrows1, rowduals_);
     1255
     1256  rcosts_ = new double[ncols0_];
     1257  ClpDisjointCopyN(si->getReducedCost(), ncols1, rcosts_);
     1258  if (maxmin_<0.0) {
     1259    // change so will look as if minimize
     1260    int i;
     1261    for (i=0;i<nrows1;i++)
     1262      rowduals_[i] = - rowduals_[i];
     1263    for (i=0;i<ncols1;i++) {
     1264      rcosts_[i] = - rcosts_[i];
     1265    }
     1266  }
     1267
     1268  //ClpDisjointCopyN(si->getRowUpper(), nrows1, rup_);
     1269  //ClpDisjointCopyN(si->getRowLower(), nrows1, rlo_);
     1270
     1271  ClpDisjointCopyN(si->getColSolution(), ncols1, sol_);
     1272  si->setDblParam(ClpObjOffset,originalOffset_);
     1273
     1274  for (int j=0; j<ncols1; j++) {
     1275    CoinBigIndex kcs = mcstrt_[j];
     1276    CoinBigIndex kce = kcs + hincol_[j];
     1277    for (CoinBigIndex k=kcs; k<kce; ++k) {
     1278      link_[k] = k+1;
     1279    }
     1280  }
     1281  {
     1282    int ml = maxlink_;
     1283    for (CoinBigIndex k=nelemsr; k<ml; ++k)
     1284      link_[k] = k+1;
     1285    link_[ml-1] = NO_LINK;
     1286  }
     1287  free_list_ = nelemsr;
     1288}
     1289/* This is main part of Presolve */
     1290ClpSimplex *
     1291ClpPresolve::gutsOfPresolvedModel(ClpSimplex * originalModel,
     1292                                  double feasibilityTolerance,
     1293                                  bool keepIntegers,
     1294                                  int numberPasses,
     1295                                  bool dropNames)
     1296{
     1297  ncols_ = originalModel->getNumCols();
     1298  nrows_ = originalModel->getNumRows();
     1299  nelems_ = originalModel->getNumElements();
    861300  numberPasses_ = numberPasses;
    871301
    88   double maxmin = si.getObjSense();
    89   originalModel_ = &si;
     1302  double maxmin = originalModel->getObjSense();
     1303  originalModel_ = originalModel;
    901304  delete [] originalColumn_;
    911305  originalColumn_ = new int[ncols_];
    921306  delete [] originalRow_;
    931307  originalRow_ = new int[nrows_];
    94 
    95   // Check matrix
    96   if (!si.clpMatrix()->allElementsInRange(&si,si.getSmallElementValue(),
    97                                           1.0e20))
    98     return NULL;
    991308
    1001309  // result is 0 - okay, 1 infeasible, -1 go round again
     
    1041313  presolvedModel_=NULL;
    1051314  // Messages
    106   CoinMessages messages = CoinMessage(si.messages().language());
     1315  CoinMessages messages = CoinMessage(originalModel->messages().language());
    1071316  while (result==-1) {
    1081317
    1091318    // make new copy
    110     delete presolvedModel_;
    111     presolvedModel_ = new ClpSimplex(si);
     1319    if (saveFile_=="") {
     1320      delete presolvedModel_;
     1321      presolvedModel_ = new ClpSimplex(*originalModel);
     1322    } else {
     1323      presolvedModel_=originalModel;
     1324    }
    1121325    presolvedModel_->dropNames();
    1131326
     
    2201433      delete [] prob.originalRow_;
    2211434      prob.originalRow_=NULL;
    222       if (!dropNames&&si.lengthNames()) {
     1435      if (!dropNames&&originalModel->lengthNames()) {
    2231436        // Redo names
    2241437        int iRow;
     
    2271440        for (iRow=0;iRow<nrowsNow;iRow++) {
    2281441          int kRow = originalRow_[iRow];
    229           rowNames.push_back(si.rowName(kRow));
     1442          rowNames.push_back(originalModel->rowName(kRow));
    2301443        }
    2311444     
     
    2351448        for (iColumn=0;iColumn<ncolsNow;iColumn++) {
    2361449          int kColumn = originalColumn_[iColumn];
    237           columnNames.push_back(si.columnName(kColumn));
     1450          columnNames.push_back(originalModel->columnName(kColumn));
    2381451        }
    2391452        presolvedModel_->copyNames(rowNames,columnNames);
     
    3141527}
    3151528
    316 // Return pointer to presolved model
    317 ClpSimplex *
    318 ClpPresolve::model() const
    319 {
    320   return presolvedModel_;
    321 }
    322 // Return pointer to original model
    323 ClpSimplex *
    324 ClpPresolve::originalModel() const
    325 {
    326   return originalModel_;
    327 }
    328 void
    329 ClpPresolve::postsolve(bool updateStatus)
    330 {
    331   // Messages
    332   CoinMessages messages = CoinMessage(presolvedModel_->messages().language());
    333   if (!presolvedModel_->isProvenOptimal()) {
    334     presolvedModel_->messageHandler()->message(COIN_PRESOLVE_NONOPTIMAL,
    335                                              messages)
    336                                                <<CoinMessageEol;
    337   }
    338 
    339   // this is the size of the original problem
    340   const int ncols0  = ncols_;
    341   const int nrows0  = nrows_;
    342   const CoinBigIndex nelems0 = nelems_;
    343 
    344   // reality check
    345   assert(ncols0==originalModel_->getNumCols());
    346   assert(nrows0==originalModel_->getNumRows());
    347 
    348   // this is the reduced problem
    349   int ncols = presolvedModel_->getNumCols();
    350   int nrows = presolvedModel_->getNumRows();
    351 
    352   double *acts = originalModel_->primalRowSolution();
    353   double *sol  = originalModel_->primalColumnSolution();
    354 
    355   unsigned char * rowstat=NULL;
    356   unsigned char * colstat = NULL;
    357   if (updateStatus) {
    358     unsigned char *status = originalModel_->statusArray();
    359     rowstat = status + ncols0;
    360     colstat = status;
    361     memcpy(colstat, presolvedModel_->statusArray(), ncols);
    362     memcpy(rowstat, presolvedModel_->statusArray()+ncols, nrows);
    363   }
    364 
    365 
    366   CoinPostsolveMatrix prob(presolvedModel_,
    367                        ncols0,
    368                        nrows0,
    369                        nelems0,
    370                        presolvedModel_->getObjSense(),
    371                        // end prepost
    372                        
    373                        sol, acts,
    374                        colstat, rowstat);
    375    
    376   postsolve(prob);
    377 
    378 }
    379 
    380 // return pointer to original columns
    381 const int *
    382 ClpPresolve::originalColumns() const
    383 {
    384   return originalColumn_;
    385 }
    386 // return pointer to original rows
    387 const int *
    388 ClpPresolve::originalRows() const
    389 {
    390   return originalRow_;
    391 }
    392 // Set pointer to original model
    393 void
    394 ClpPresolve::setOriginalModel(ClpSimplex * model)
    395 {
    396   originalModel_=model;
    397 }
    398 #if 0
    399 // A lazy way to restrict which transformations are applied
    400 // during debugging.
    401 static int ATOI(const char *name)
    402 {
    403  return true;
    404 #if     DEBUG_PRESOLVE || PRESOLVE_SUMMARY
    405   if (getenv(name)) {
    406     int val = atoi(getenv(name));
    407     printf("%s = %d\n", name, val);
    408     return (val);
    409   } else {
    410     if (strcmp(name,"off"))
    411       return (true);
    412     else
    413       return (false);
    414   }
    415 #else
    416   return (true);
    417 #endif
    418 }
    419 #endif
    420 //#define DEBUG_PRESOLVE 1
    421 #if DEBUG_PRESOLVE
    422 void check_sol(CoinPresolveMatrix *prob,double tol)
    423 {
    424   double *colels        = prob->colels_;
    425   int *hrow             = prob->hrow_;
    426   int *mcstrt           = prob->mcstrt_;
    427   int *hincol           = prob->hincol_;
    428   int *hinrow           = prob->hinrow_;
    429   int ncols             = prob->ncols_;
    430 
    431 
    432   double * csol = prob->sol_;
    433   double * acts = prob->acts_;
    434   double * clo = prob->clo_;
    435   double * cup = prob->cup_;
    436   int nrows = prob->nrows_;
    437   double * rlo = prob->rlo_;
    438   double * rup = prob->rup_;
    439 
    440   int colx;
    441 
    442   double * rsol = new double[nrows];
    443   memset(rsol,0,nrows*sizeof(double));
    444 
    445   for (colx = 0; colx < ncols; ++colx) {
    446     if (1) {
    447       CoinBigIndex k = mcstrt[colx];
    448       int nx = hincol[colx];
    449       double solutionValue = csol[colx];
    450       for (int i=0; i<nx; ++i) {
    451         int row = hrow[k];
    452         double coeff = colels[k];
    453         k++;
    454         rsol[row] += solutionValue*coeff;
    455       }
    456       if (csol[colx]<clo[colx]-tol) {
    457         printf("low CSOL:  %d  - %g %g %g\n",
    458                    colx, clo[colx], csol[colx], cup[colx]);
    459       } else if (csol[colx]>cup[colx]+tol) {
    460         printf("high CSOL:  %d  - %g %g %g\n",
    461                    colx, clo[colx], csol[colx], cup[colx]);
    462       }
    463     }
    464   }
    465   int rowx;
    466   for (rowx = 0; rowx < nrows; ++rowx) {
    467     if (hinrow[rowx]) {
    468       if (fabs(rsol[rowx]-acts[rowx])>tol)
    469         printf("inacc RSOL:  %d - %g %g (acts_ %g) %g\n",
    470                    rowx,  rlo[rowx], rsol[rowx], acts[rowx], rup[rowx]);
    471       if (rsol[rowx]<rlo[rowx]-tol) {
    472         printf("low RSOL:  %d - %g %g %g\n",
    473                    rowx,  rlo[rowx], rsol[rowx], rup[rowx]);
    474       } else if (rsol[rowx]>rup[rowx]+tol ) {
    475         printf("high RSOL:  %d - %g %g %g\n",
    476                    rowx,  rlo[rowx], rsol[rowx], rup[rowx]);
    477       }
    478     }
    479   }
    480   delete [] rsol;
    481 }
    482 #endif
    483 // This is the presolve loop.
    484 // It is a separate virtual function so that it can be easily
    485 // customized by subclassing CoinPresolve.
    486 const CoinPresolveAction *ClpPresolve::presolve(CoinPresolveMatrix *prob)
    487 {
    488   paction_ = 0;
    489 
    490   prob->status_=0; // say feasible
    491 
    492   paction_ = make_fixed(prob, paction_);
    493   // if integers then switch off dual stuff
    494   // later just do individually
    495   bool doDualStuff = (presolvedModel_->integerInformation()==NULL);
    496 
    497 #if     CHECK_CONSISTENCY
    498   presolve_links_ok(prob->rlink_, prob->mrstrt_, prob->hinrow_, prob->nrows_);
    499 #endif
    500 
    501   if (!prob->status_) {
    502 #if 0
    503     const bool slackd = ATOI("SLACKD")!=0;
    504     //const bool forcing = ATOI("FORCING")!=0;
    505     const bool doubleton = ATOI("DOUBLETON")!=0;
    506     const bool forcing = ATOI("off")!=0;
    507     const bool ifree = ATOI("off")!=0;
    508     const bool zerocost = ATOI("off")!=0;
    509     const bool dupcol = ATOI("off")!=0;
    510     const bool duprow = ATOI("off")!=0;
    511     const bool dual = ATOI("off")!=0;
    512 #else
    513     // normal
    514 #if 0
    515     const bool slackd = true;
    516     const bool doubleton = false;
    517     const bool tripleton = true;
    518     const bool forcing = false;
    519     const bool ifree = true;
    520     const bool zerocost = true;
    521     const bool dupcol = false;
    522     const bool duprow = false;
    523     const bool dual = doDualStuff;
    524 #else
    525     const bool slackd = true;
    526     const bool doubleton = true;
    527     const bool tripleton = true;
    528     const bool forcing = true;
    529     const bool ifree = true;
    530     const bool zerocost = true;
    531     const bool dupcol = true;
    532     const bool duprow = true;
    533     const bool dual = doDualStuff;
    534 #endif
    535 #endif
    536    
    537     // some things are expensive so just do once (normally)
    538 
    539     int i;
    540     // say look at all
    541     if (!prob->anyProhibited()) {
    542       for (i=0;i<nrows_;i++)
    543         prob->rowsToDo_[i]=i;
    544       prob->numberRowsToDo_=nrows_;
    545       for (i=0;i<ncols_;i++)
    546         prob->colsToDo_[i]=i;
    547       prob->numberColsToDo_=ncols_;
    548     } else {
    549       // some stuff must be left alone
    550       prob->numberRowsToDo_=0;
    551       for (i=0;i<nrows_;i++)
    552         if (!prob->rowProhibited(i))
    553             prob->rowsToDo_[prob->numberRowsToDo_++]=i;
    554       prob->numberColsToDo_=0;
    555       for (i=0;i<ncols_;i++)
    556         if (!prob->colProhibited(i))
    557             prob->colsToDo_[prob->numberColsToDo_++]=i;
    558     }
    559 
    560 
    561     int iLoop;
    562 #if     DEBUG_PRESOLVE
    563     check_sol(prob,1.0e0);
    564 #endif
    565 
    566     // Check number rows dropped
    567     int lastDropped=0;
    568     for (iLoop=0;iLoop<numberPasses_;iLoop++) {
    569 #ifdef PRESOLVE_SUMMARY
    570       printf("Starting major pass %d\n",iLoop+1);
    571 #endif
    572       const CoinPresolveAction * const paction0 = paction_;
    573       // look for substitutions with no fill
    574       int fill_level=2;
    575       //fill_level=10;
    576       //printf("** fill_level == 10 !\n");
    577       int whichPass=0;
    578       while (1) {
    579         whichPass++;
    580         const CoinPresolveAction * const paction1 = paction_;
    581 
    582         if (slackd) {
    583           bool notFinished = true;
    584           while (notFinished)
    585             paction_ = slack_doubleton_action::presolve(prob, paction_,
    586                                                         notFinished);
    587           if (prob->status_)
    588             break;
    589         }
    590 
    591         if (doubleton) {
    592           paction_ = doubleton_action::presolve(prob, paction_);
    593           if (prob->status_)
    594             break;
    595         }
    596 
    597         if (tripleton) {
    598           paction_ = tripleton_action::presolve(prob, paction_);
    599           if (prob->status_)
    600             break;
    601         }
    602 
    603         if (zerocost) {
    604           paction_ = do_tighten_action::presolve(prob, paction_);
    605           if (prob->status_)
    606             break;
    607         }
    608 
    609         if (forcing) {
    610           paction_ = forcing_constraint_action::presolve(prob, paction_);
    611           if (prob->status_)
    612             break;
    613         }
    614 
    615         if (ifree) {
    616           paction_ = implied_free_action::presolve(prob, paction_,fill_level);
    617           if (prob->status_)
    618             break;
    619         }
    620 
    621 #if     DEBUG_PRESOLVE
    622         check_sol(prob,1.0e0);
    623 #endif
    624 
    625 #if     CHECK_CONSISTENCY
    626         presolve_links_ok(prob->rlink_, prob->mrstrt_, prob->hinrow_,
    627                           prob->nrows_);
    628 #endif
    629 
    630 #if     DEBUG_PRESOLVE
    631         presolve_no_zeros(prob->mcstrt_, prob->colels_, prob->hincol_,
    632                           prob->ncols_);
    633 #endif
    634 #if     CHECK_CONSISTENCY
    635         prob->consistent();
    636 #endif
    637 
    638          
    639         // set up for next pass
    640         // later do faster if many changes i.e. memset and memcpy
    641         prob->numberRowsToDo_ = prob->numberNextRowsToDo_;
    642         int kcheck;
    643         bool found=false;
    644         kcheck=-1;
    645         for (i=0;i<prob->numberNextRowsToDo_;i++) {
    646           int index = prob->nextRowsToDo_[i];
    647           prob->unsetRowChanged(index);
    648           prob->rowsToDo_[i] = index;
    649           if (index==kcheck) {
    650             printf("row %d on list after pass %d\n",kcheck,
    651                    whichPass);
    652             found=true;
    653           }
    654         }
    655         if (!found&&kcheck>=0)
    656           prob->rowsToDo_[prob->numberRowsToDo_++]=kcheck;
    657         prob->numberNextRowsToDo_=0;
    658         prob->numberColsToDo_ = prob->numberNextColsToDo_;
    659         kcheck=-1;
    660         found=false;
    661         for (i=0;i<prob->numberNextColsToDo_;i++) {
    662           int index = prob->nextColsToDo_[i];
    663           prob->unsetColChanged(index);
    664           prob->colsToDo_[i] = index;
    665           if (index==kcheck) {
    666             printf("col %d on list after pass %d\n",kcheck,
    667                    whichPass);
    668             found=true;
    669           }
    670         }
    671         if (!found&&kcheck>=0)
    672           prob->colsToDo_[prob->numberColsToDo_++]=kcheck;
    673         prob->numberNextColsToDo_=0;
    674         if (paction_ == paction1&&fill_level>0)
    675           break;
    676       }
    677       // say look at all
    678       int i;
    679       if (!prob->anyProhibited()) {
    680         for (i=0;i<nrows_;i++)
    681           prob->rowsToDo_[i]=i;
    682         prob->numberRowsToDo_=nrows_;
    683         for (i=0;i<ncols_;i++)
    684           prob->colsToDo_[i]=i;
    685         prob->numberColsToDo_=ncols_;
    686       } else {
    687         // some stuff must be left alone
    688         prob->numberRowsToDo_=0;
    689         for (i=0;i<nrows_;i++)
    690           if (!prob->rowProhibited(i))
    691             prob->rowsToDo_[prob->numberRowsToDo_++]=i;
    692         prob->numberColsToDo_=0;
    693         for (i=0;i<ncols_;i++)
    694           if (!prob->colProhibited(i))
    695             prob->colsToDo_[prob->numberColsToDo_++]=i;
    696       }
    697       // now expensive things
    698       // this caused world.mps to run into numerical difficulties
    699 #ifdef PRESOLVE_SUMMARY
    700       printf("Starting expensive\n");
    701 #endif
    702 
    703       if (dual) {
    704         int itry;
    705         for (itry=0;itry<5;itry++) {
    706           const CoinPresolveAction * const paction2 = paction_;
    707           paction_ = remove_dual_action::presolve(prob, paction_);
    708           if (prob->status_)
    709             break;
    710           if (ifree) {
    711             int fill_level=0; // switches off substitution
    712             paction_ = implied_free_action::presolve(prob, paction_,fill_level);
    713             if (prob->status_)
    714               break;
    715           }
    716           if (paction_ == paction2)
    717             break;
    718         }
    719       }
    720 #if     DEBUG_PRESOLVE
    721       check_sol(prob,1.0e0);
    722 #endif
    723       if (dupcol) {
    724         paction_ = dupcol_action::presolve(prob, paction_);
    725         if (prob->status_)
    726           break;
    727       }
    728 #if     DEBUG_PRESOLVE
    729         check_sol(prob,1.0e0);
    730 #endif
    731      
    732       if (duprow) {
    733         paction_ = duprow_action::presolve(prob, paction_);
    734         if (prob->status_)
    735           break;
    736       }
    737 #if     DEBUG_PRESOLVE
    738       check_sol(prob,1.0e0);
    739 #endif
    740       {
    741         int * hinrow = prob->hinrow_;
    742         int numberDropped=0;
    743         for (i=0;i<nrows_;i++)
    744           if (!hinrow[i])
    745             numberDropped++;
    746         //printf("%d rows dropped after pass %d\n",numberDropped,
    747         //     iLoop+1);
    748         if (numberDropped==lastDropped)
    749           break;
    750         else
    751           lastDropped = numberDropped;
    752       }
    753       if (paction_ == paction0)
    754         break;
    755          
    756     }
    757   }
    758   if (!prob->status_) {
    759     paction_ = drop_zero_coefficients(prob, paction_);
    760 #if     DEBUG_PRESOLVE
    761         check_sol(prob,1.0e0);
    762 #endif
    763 
    764     paction_ = drop_empty_cols_action::presolve(prob, paction_);
    765     paction_ = drop_empty_rows_action::presolve(prob, paction_);
    766 #if     DEBUG_PRESOLVE
    767         check_sol(prob,1.0e0);
    768 #endif
    769   }
    770  
    771   // Messages
    772   CoinMessages messages = CoinMessage(prob->messages().language());
    773   if (prob->status_) {
    774     if (prob->status_==1)
    775           prob->messageHandler()->message(COIN_PRESOLVE_INFEAS,
    776                                              messages)
    777                                                <<prob->feasibilityTolerance_
    778                                                <<CoinMessageEol;
    779     else if (prob->status_==2)
    780           prob->messageHandler()->message(COIN_PRESOLVE_UNBOUND,
    781                                              messages)
    782                                                <<CoinMessageEol;
    783     else
    784           prob->messageHandler()->message(COIN_PRESOLVE_INFEASUNBOUND,
    785                                              messages)
    786                                                <<CoinMessageEol;
    787     // get rid of data
    788     gutsOfDestroy();
    789   }
    790   return (paction_);
    791 }
    792 
    793 void check_djs(CoinPostsolveMatrix *prob);
    794 
    795 
    796 // We could have implemented this by having each postsolve routine
    797 // directly call the next one, but this may make it easier to add debugging checks.
    798 void ClpPresolve::postsolve(CoinPostsolveMatrix &prob)
    799 {
    800   const CoinPresolveAction *paction = paction_;
    801 
    802   if (prob.colstat_)
    803     prob.check_nbasic();
    804  
    805 #if     DEBUG_PRESOLVE
    806   check_djs(&prob);
    807 #endif
    808  
    809  
    810   while (paction) {
    811 #if     DEBUG_PRESOLVE
    812     printf("POSTSOLVING %s\n", paction->name());
    813 #endif
    814 
    815     paction->postsolve(&prob);
    816    
    817 #if     DEBUG_PRESOLVE
    818     if (prob.colstat_)
    819       prob.check_nbasic();
    820 #endif
    821     paction = paction->next;
    822 #if     DEBUG_PRESOLVE
    823     check_djs(&prob);
    824 #endif
    825   }   
    826  
    827 #if     0 && DEBUG_PRESOLVE
    828   for (i=0; i<ncols0; i++) {
    829     if (!cdone[i]) {
    830       printf("!cdone[%d]\n", i);
    831       abort();
    832     }
    833   }
    834  
    835   for (i=0; i<nrows0; i++) {
    836     if (!rdone[i]) {
    837       printf("!rdone[%d]\n", i);
    838       abort();
    839     }
    840   }
    841  
    842  
    843   for (i=0; i<ncols0; i++) {
    844     if (sol[i] < -1e10 || sol[i] > 1e10)
    845       printf("!!!%d %g\n", i, sol[i]);
    846    
    847   }
    848  
    849  
    850 #endif
    851  
    852 #if     0 && DEBUG_PRESOLVE
    853   // debug check:  make sure we ended up with same original matrix
    854   {
    855     int identical = 1;
    856    
    857     for (int i=0; i<ncols0; i++) {
    858       PRESOLVEASSERT(hincol[i] == &prob->mcstrt0[i+1] - &prob->mcstrt0[i]);
    859       CoinBigIndex kcs0 = &prob->mcstrt0[i];
    860       CoinBigIndex kcs = mcstrt[i];
    861       int n = hincol[i];
    862       for (int k=0; k<n; k++) {
    863         CoinBigIndex k1 = presolve_find_row1(&prob->hrow0[kcs0+k], kcs, kcs+n, hrow);
    864 
    865         if (k1 == kcs+n) {
    866           printf("ROW %d NOT IN COL %d\n", &prob->hrow0[kcs0+k], i);
    867           abort();
    868         }
    869 
    870         if (colels[k1] != &prob->dels0[kcs0+k])
    871           printf("BAD COLEL[%d %d %d]:  %g\n",
    872                  k1, i, &prob->hrow0[kcs0+k], colels[k1] - &prob->dels0[kcs0+k]);
    873 
    874         if (kcs0+k != k1)
    875           identical=0;
    876       }
    877     }
    878     printf("identical? %d\n", identical);
    879   }
    880 #endif
    881   // put back duals
    882   memcpy(originalModel_->dualRowSolution(),prob.rowduals_,
    883          nrows_*sizeof(double));
    884   double maxmin = originalModel_->getObjSense();
    885   if (maxmin<0.0) {
    886     // swap signs
    887     int i;
    888     double * pi = originalModel_->dualRowSolution();
    889     for (i=0;i<nrows_;i++)
    890       pi[i] = -pi[i];
    891   }
    892   // Now check solution
    893   memcpy(originalModel_->dualColumnSolution(),
    894          originalModel_->objective(),ncols_*sizeof(double));
    895   originalModel_->transposeTimes(-1.0,
    896                                  originalModel_->dualRowSolution(),
    897                                  originalModel_->dualColumnSolution());
    898   memset(originalModel_->primalRowSolution(),0,nrows_*sizeof(double));
    899   originalModel_->times(1.0,originalModel_->primalColumnSolution(),
    900                         originalModel_->primalRowSolution());
    901   originalModel_->checkSolution();
    902   // Messages
    903   CoinMessages messages = CoinMessage(presolvedModel_->messages().language());
    904   presolvedModel_->messageHandler()->message(COIN_PRESOLVE_POSTSOLVE,
    905                                             messages)
    906                                               <<originalModel_->objectiveValue()
    907                                               <<originalModel_->sumDualInfeasibilities()
    908                                               <<originalModel_->numberDualInfeasibilities()
    909                                               <<originalModel_->sumPrimalInfeasibilities()
    910                                               <<originalModel_->numberPrimalInfeasibilities()
    911                                                <<CoinMessageEol;
    912  
    913   //originalModel_->objectiveValue_=objectiveValue_;
    914   originalModel_->setNumberIterations(presolvedModel_->numberIterations());
    915   if (!presolvedModel_->status()) {
    916     if (!originalModel_->numberDualInfeasibilities()&&
    917         !originalModel_->numberPrimalInfeasibilities()) {
    918       originalModel_->setProblemStatus( 0);
    919     } else {
    920       originalModel_->setProblemStatus( -1);
    921       presolvedModel_->messageHandler()->message(COIN_PRESOLVE_NEEDS_CLEANING,
    922                                             messages)
    923                                               <<CoinMessageEol;
    924     }
    925   } else {
    926     originalModel_->setProblemStatus( presolvedModel_->status());
    927   }
    928 }
    929 
    930 
    931 
    932 
    933 
    934 
    935 
    936 
    937 #if     DEBUG_PRESOLVE
    938 void check_djs(CoinPostsolveMatrix *prob)
    939 {
    940   //return;
    941   double *colels        = prob->colels_;
    942   int *hrow             = prob->hrow_;
    943   int *mcstrt           = prob->mcstrt_;
    944   int *hincol           = prob->hincol_;
    945   int *link             = prob->link_;
    946   int ncols             = prob->ncols_;
    947 
    948   double *dcost = prob->cost_;
    949 
    950   double *rcosts        = prob->rcosts_;
    951 
    952   double *rowduals = prob->rowduals_;
    953 
    954   const double maxmin   = prob->maxmin_;
    955 
    956   char *cdone   = prob->cdone_;
    957 
    958   double * csol = prob->sol_;
    959   double * clo = prob->clo_;
    960   double * cup = prob->cup_;
    961   int nrows = prob->nrows_;
    962   double * rlo = prob->rlo_;
    963   double * rup = prob->rup_;
    964   char *rdone   = prob->rdone_;
    965 
    966   int colx;
    967 
    968   double * rsol = new double[nrows];
    969   memset(rsol,0,nrows*sizeof(double));
    970 
    971   for (colx = 0; colx < ncols; ++colx) {
    972     if (cdone[colx]) {
    973       CoinBigIndex k = mcstrt[colx];
    974       int nx = hincol[colx];
    975       double dj = maxmin * dcost[colx];
    976       double solutionValue = csol[colx];
    977       for (int i=0; i<nx; ++i) {
    978         int row = hrow[k];
    979         double coeff = colels[k];
    980         k = link[k];
    981         dj -= rowduals[row] * coeff;
    982         rsol[row] += solutionValue*coeff;
    983       }
    984       if (! (fabs(rcosts[colx] - dj) < 1.0e-4))
    985         printf("BAD DJ:  %d %g %g\n",
    986                colx, rcosts[colx], dj);
    987       if (cup[colx]-clo[colx]>1.0e-6) {
    988         if (csol[colx]<clo[colx]+1.0e-6) {
    989           if (dj <-1.0e-6)
    990             printf("neg DJ:  %d %g  - %g %g %g\n",
    991                    colx, dj, clo[colx], csol[colx], cup[colx]);
    992         } else if (csol[colx]>cup[colx]-1.0e-6) {
    993           if (dj > 1.0e-6)
    994             printf("pos DJ:  %d %g  - %g %g %g\n",
    995                    colx, dj, clo[colx], csol[colx], cup[colx]);
    996         } else {
    997           if (fabs(dj) >1.0e-6)
    998             printf("nonzero DJ:  %d %g  - %g %g %g\n",
    999                    colx, dj, clo[colx], csol[colx], cup[colx]);
    1000         }
    1001       }
    1002     }
    1003   }
    1004   int rowx;
    1005   for (rowx = 0; rowx < nrows; ++rowx) {
    1006     if (rdone[rowx]) {
    1007       if (rup[rowx]-rlo[rowx]>1.0e-6) {
    1008         double dj = rowduals[rowx];
    1009         if (rsol[rowx]<rlo[rowx]+1.0e-6) {
    1010           if (dj <-1.0e-5)
    1011             printf("neg rDJ:  %d %g  - %g %g %g\n",
    1012                    rowx, dj, rlo[rowx], rsol[rowx], rup[rowx]);
    1013         } else if (rsol[rowx]>rup[rowx]-1.0e-6) {
    1014           if (dj > 1.0e-5)
    1015             printf("pos rDJ:  %d %g  - %g %g %g\n",
    1016                    rowx, dj, rlo[rowx], rsol[rowx], rup[rowx]);
    1017         } else {
    1018           if (fabs(dj) >1.0e-5)
    1019             printf("nonzero rDJ:  %d %g  - %g %g %g\n",
    1020                    rowx, dj, rlo[rowx], rsol[rowx], rup[rowx]);
    1021         }
    1022       }
    1023     }
    1024   }
    1025   delete [] rsol;
    1026 }
    1027 #endif
    1028 
    1029 static inline double getTolerance(const ClpSimplex  *si, ClpDblParam key)
    1030 {
    1031   double tol;
    1032   if (! si->getDblParam(key, tol)) {
    1033     CoinPresolveAction::throwCoinError("getDblParam failed",
    1034                                       "CoinPrePostsolveMatrix::CoinPrePostsolveMatrix");
    1035   }
    1036   return (tol);
    1037 }
    1038 
    1039 
    1040 // Assumptions:
    1041 // 1. nrows>=m.getNumRows()
    1042 // 2. ncols>=m.getNumCols()
    1043 //
    1044 // In presolve, these values are equal.
    1045 // In postsolve, they may be inequal, since the reduced problem
    1046 // may be smaller, but we need room for the large problem.
    1047 // ncols may be larger than si.getNumCols() in postsolve,
    1048 // this at that point si will be the reduced problem,
    1049 // but we need to reserve enough space for the original problem.
    1050 CoinPrePostsolveMatrix::CoinPrePostsolveMatrix(const ClpSimplex * si,
    1051                                              int ncols_in,
    1052                                              int nrows_in,
    1053                                              CoinBigIndex nelems_in) :
    1054   ncols_(si->getNumCols()),
    1055   ncols0_(ncols_in),
    1056   nelems_(si->getNumElements()),
    1057 
    1058   mcstrt_(new CoinBigIndex[ncols_in+1]),
    1059   hincol_(new int[ncols_in+1]),
    1060   hrow_  (new int   [2*nelems_in]),
    1061   colels_(new double[2*nelems_in]),
    1062 
    1063   cost_(new double[ncols_in]),
    1064   clo_(new double[ncols_in]),
    1065   cup_(new double[ncols_in]),
    1066   rlo_(new double[nrows_in]),
    1067   rup_(new double[nrows_in]),
    1068   originalColumn_(new int[ncols_in]),
    1069   originalRow_(new int[nrows_in]),
    1070 
    1071   ztolzb_(getTolerance(si, ClpPrimalTolerance)),
    1072   ztoldj_(getTolerance(si, ClpDualTolerance)),
    1073 
    1074   maxmin_(si->getObjSense())
    1075 
    1076 {
    1077   si->getDblParam(ClpObjOffset,originalOffset_);
    1078   int ncols = si->getNumCols();
    1079   int nrows = si->getNumRows();
    1080 
    1081   ClpDisjointCopyN(si->getColLower(), ncols, clo_);
    1082   ClpDisjointCopyN(si->getColUpper(), ncols, cup_);
    1083   ClpDisjointCopyN(si->getObjCoefficients(), ncols, cost_);
    1084   ClpDisjointCopyN(si->getRowLower(), nrows,  rlo_);
    1085   ClpDisjointCopyN(si->getRowUpper(), nrows,  rup_);
    1086   int i;
    1087   for (i=0;i<ncols_in;i++)
    1088     originalColumn_[i]=i;
    1089   for (i=0;i<nrows_in;i++)
    1090     originalRow_[i]=i;
    1091   sol_=NULL;
    1092   rowduals_=NULL;
    1093   acts_=NULL;
    1094 
    1095   rcosts_=NULL;
    1096   colstat_=NULL;
    1097   rowstat_=NULL;
    1098 }
    1099 
    1100 // I am not familiar enough with CoinPackedMatrix to be confident
    1101 // that I will implement a row-ordered version of toColumnOrderedGapFree
    1102 // properly.
    1103 static bool isGapFree(const CoinPackedMatrix& matrix)
    1104 {
    1105   const CoinBigIndex * start = matrix.getVectorStarts();
    1106   const int * length = matrix.getVectorLengths();
    1107   int i;
    1108   for (i = matrix.getSizeVectorLengths() - 1; i >= 0; --i) {
    1109     if (start[i+1] - start[i] != length[i])
    1110       break;
    1111   }
    1112   return (! (i >= 0));
    1113 }
    1114 #if     DEBUG_PRESOLVE
    1115 static void matrix_bounds_ok(const double *lo, const double *up, int n)
    1116 {
    1117   int i;
    1118   for (i=0; i<n; i++) {
    1119     PRESOLVEASSERT(lo[i] <= up[i]);
    1120     PRESOLVEASSERT(lo[i] < PRESOLVE_INF);
    1121     PRESOLVEASSERT(-PRESOLVE_INF < up[i]);
    1122   }
    1123 }
    1124 #endif
    1125 CoinPresolveMatrix::CoinPresolveMatrix(int ncols0_in,
    1126                                      double maxmin_,
    1127                                      // end prepost members
    1128 
    1129                                      ClpSimplex * si,
    1130 
    1131                                      // rowrep
    1132                                      int nrows_in,
    1133                                      CoinBigIndex nelems_in,
    1134                                bool doStatus,
    1135                                double nonLinearValue) :
    1136 
    1137   CoinPrePostsolveMatrix(si,
    1138                         ncols0_in, nrows_in, nelems_in),
    1139   clink_(new presolvehlink[ncols0_in+1]),
    1140   rlink_(new presolvehlink[nrows_in+1]),
    1141 
    1142   dobias_(0.0),
    1143 
    1144   nrows_(si->getNumRows()),
    1145 
    1146   // temporary init
    1147   mrstrt_(new CoinBigIndex[nrows_in+1]),
    1148   hinrow_(new int[nrows_in+1]),
    1149   rowels_(new double[2*nelems_in]),
    1150   hcol_(new int[2*nelems_in]),
    1151   integerType_(new char[ncols0_in]),
    1152   feasibilityTolerance_(0.0),
    1153   status_(-1),
    1154   rowsToDo_(new int [nrows_in]),
    1155   numberRowsToDo_(0),
    1156   nextRowsToDo_(new int[nrows_in]),
    1157   numberNextRowsToDo_(0),
    1158   colsToDo_(new int [ncols0_in]),
    1159   numberColsToDo_(0),
    1160   nextColsToDo_(new int[ncols0_in]),
    1161   numberNextColsToDo_(0)
    1162 
    1163 {
    1164   const int bufsize = 2*nelems_in;
    1165 
    1166   // Set up change bits etc
    1167   rowChanged_ = new unsigned char[nrows_];
    1168   memset(rowChanged_,0,nrows_);
    1169   colChanged_ = new unsigned char[ncols_];
    1170   memset(colChanged_,0,ncols_);
    1171   CoinPackedMatrix * m = si->matrix();
    1172 
    1173   // The coefficient matrix is a big hunk of stuff.
    1174   // Do the copy here to try to avoid running out of memory.
    1175 
    1176   const CoinBigIndex * start = m->getVectorStarts();
    1177   const int * length = m->getVectorLengths();
    1178   const int * row = m->getIndices();
    1179   const double * element = m->getElements();
    1180   int icol,nel=0;
    1181   mcstrt_[0]=0;
    1182   for (icol=0;icol<ncols_;icol++) {
    1183     int j;
    1184     for (j=start[icol];j<start[icol]+length[icol];j++) {
    1185       hrow_[nel]=row[j];
    1186       colels_[nel++]=element[j];
    1187     }
    1188     mcstrt_[icol+1]=nel;
    1189   }
    1190   assert(mcstrt_[ncols_] == nelems_);
    1191   ClpDisjointCopyN(m->getVectorLengths(),ncols_,  hincol_);
    1192 
    1193   // same thing for row rep
    1194   m = new CoinPackedMatrix();
    1195   m->reverseOrderedCopyOf(*si->matrix());
    1196   m->removeGaps();
    1197 
    1198 
    1199   ClpDisjointCopyN(m->getVectorStarts(),  nrows_,  mrstrt_);
    1200   mrstrt_[nrows_] = nelems_;
    1201   ClpDisjointCopyN(m->getVectorLengths(), nrows_,  hinrow_);
    1202   ClpDisjointCopyN(m->getIndices(),       nelems_, hcol_);
    1203   ClpDisjointCopyN(m->getElements(),      nelems_, rowels_);
    1204 
    1205   delete m;
    1206   if (si->integerInformation()) {
    1207     memcpy(integerType_,si->integerInformation(),ncols_*sizeof(char));
    1208   } else {
    1209     ClpFillN<char>(integerType_, ncols_, 0);
    1210   }
    1211 
    1212   // Set up prohibited bits if needed
    1213   if (nonLinearValue) {
    1214     anyProhibited_ = true;
    1215     for (icol=0;icol<ncols_;icol++) {
    1216       int j;
    1217       bool nonLinearColumn = false;
    1218       if (cost_[icol]==nonLinearValue)
    1219         nonLinearColumn=true;
    1220       for (j=mcstrt_[icol];j<mcstrt_[icol+1];j++) {
    1221         if (colels_[j]==nonLinearValue) {
    1222           nonLinearColumn=true;
    1223           setRowProhibited(hrow_[j]);
    1224         }
    1225       }
    1226       if (nonLinearColumn)
    1227         setColProhibited(icol);
    1228     }
    1229   } else {
    1230     anyProhibited_ = false;
    1231   }
    1232 
    1233   if (doStatus) {
    1234     // allow for status and solution
    1235     sol_ = new double[ncols_];
    1236     memcpy(sol_,si->primalColumnSolution(),ncols_*sizeof(double));;
    1237     acts_ = new double [nrows_];
    1238     memcpy(acts_,si->primalRowSolution(),nrows_*sizeof(double));
    1239     if (!si->statusArray())
    1240       si->createStatus();
    1241     colstat_ = new unsigned char [nrows_+ncols_];
    1242     memcpy(colstat_,si->statusArray(),
    1243            (nrows_+ncols_)*sizeof(unsigned char));
    1244     rowstat_ = colstat_+ncols_;
    1245   }
    1246 
    1247   // the original model's fields are now unneeded - free them
    1248  
    1249   si->resize(0,0);
    1250 
    1251 #if     DEBUG_PRESOLVE
    1252   matrix_bounds_ok(rlo_, rup_, nrows_);
    1253   matrix_bounds_ok(clo_, cup_, ncols_);
    1254 #endif
    1255 
    1256 #if 0
    1257   for (i=0; i<nrows; ++i)
    1258     printf("NR: %6d\n", hinrow[i]);
    1259   for (int i=0; i<ncols; ++i)
    1260     printf("NC: %6d\n", hincol[i]);
    1261 #endif
    1262 
    1263   presolve_make_memlists(mcstrt_, hincol_, clink_, ncols_);
    1264   presolve_make_memlists(mrstrt_, hinrow_, rlink_, nrows_);
    1265 
    1266   // this allows last col/row to expand up to bufsize-1 (22);
    1267   // this must come after the calls to presolve_prefix
    1268   mcstrt_[ncols_] = bufsize-1;
    1269   mrstrt_[nrows_] = bufsize-1;
    1270 
    1271 #if     CHECK_CONSISTENCY
    1272   consistent(false);
    1273 #endif
    1274 }
    1275 
    1276 void CoinPresolveMatrix::update_model(ClpSimplex * si,
    1277                                      int nrows0,
    1278                                      int ncols0,
    1279                                      CoinBigIndex nelems0)
    1280 {
    1281   si->loadProblem(ncols_, nrows_, mcstrt_, hrow_, colels_, hincol_,
    1282                  clo_, cup_, cost_, rlo_, rup_);
    1283 
    1284   delete [] si->integerInformation();
    1285   int numberIntegers=0;
    1286   for (int i=0; i<ncols_; i++) {
    1287     if (integerType_[i])
    1288       numberIntegers++;
    1289   }
    1290   if (numberIntegers)
    1291     si->copyInIntegerInformation(integerType_);
    1292   else
    1293     si->copyInIntegerInformation(NULL);
    1294 
    1295 #if     PRESOLVE_SUMMARY
    1296   printf("NEW NCOL/NROW/NELS:  %d(-%d) %d(-%d) %d(-%d)\n",
    1297          ncols_, ncols0-ncols_,
    1298          nrows_, nrows0-nrows_,
    1299          si->getNumElements(), nelems0-si->getNumElements());
    1300 #endif
    1301   si->setDblParam(ClpObjOffset,originalOffset_-dobias_);
    1302 
    1303 }
    1304 
    1305 
    1306 
    1307 
    1308 
    1309 
    1310 
    1311 
    1312 
    1313 
    1314 
    1315 ////////////////  POSTSOLVE
    1316 
    1317 CoinPostsolveMatrix::CoinPostsolveMatrix(ClpSimplex*  si,
    1318                                        int ncols0_in,
    1319                                        int nrows0_in,
    1320                                        CoinBigIndex nelems0,
    1321                                    
    1322                                        double maxmin_,
    1323                                        // end prepost members
    1324 
    1325                                        double *sol_in,
    1326                                        double *acts_in,
    1327 
    1328                                        unsigned char *colstat_in,
    1329                                        unsigned char *rowstat_in) :
    1330   CoinPrePostsolveMatrix(si,
    1331                         ncols0_in, nrows0_in, nelems0),
    1332 
    1333   free_list_(0),
    1334   // link, free_list, maxlink
    1335   maxlink_(2*nelems0),
    1336   link_(new int[/*maxlink*/ 2*nelems0]),
    1337      
    1338   cdone_(new char[ncols0_]),
    1339   rdone_(new char[nrows0_in]),
    1340 
    1341   nrows_(si->getNumRows()),
    1342   nrows0_(nrows0_in)
    1343 {
    1344 
    1345   sol_=sol_in;
    1346   rowduals_=NULL;
    1347   acts_=acts_in;
    1348 
    1349   rcosts_=NULL;
    1350   colstat_=colstat_in;
    1351   rowstat_=rowstat_in;
    1352 
    1353   // this is the *reduced* model, which is probably smaller
    1354   int ncols1 = si->getNumCols();
    1355   int nrows1 = si->getNumRows();
    1356 
    1357   const CoinPackedMatrix * m = si->matrix();
    1358 
    1359   if (! isGapFree(*m)) {
    1360     CoinPresolveAction::throwCoinError("Matrix not gap free",
    1361                                       "CoinPostsolveMatrix");
    1362   }
    1363 
    1364   const CoinBigIndex nelemsr = m->getNumElements();
    1365 
    1366   ClpDisjointCopyN(m->getVectorStarts(), ncols1, mcstrt_);
    1367   mcstrt_[ncols_] = nelems0;    // ??
    1368   ClpDisjointCopyN(m->getVectorLengths(),ncols1,  hincol_);
    1369   ClpDisjointCopyN(m->getIndices(),      nelemsr, hrow_);
    1370   ClpDisjointCopyN(m->getElements(),     nelemsr, colels_);
    1371 
    1372 
    1373 #if     0 && DEBUG_PRESOLVE
    1374   presolve_check_costs(model, &colcopy);
    1375 #endif
    1376 
    1377   // This determines the size of the data structure that contains
    1378   // the matrix being postsolved.  Links are taken from the free_list
    1379   // to recreate matrix entries that were presolved away,
    1380   // and links are added to the free_list when entries created during
    1381   // presolve are discarded.  There is never a need to gc this list.
    1382   // Naturally, it should contain
    1383   // exactly nelems0 entries "in use" when postsolving is done,
    1384   // but I don't know whether the matrix could temporarily get
    1385   // larger during postsolving.  Substitution into more than two
    1386   // rows could do that, in principle.  I am being very conservative
    1387   // here by reserving much more than the amount of space I probably need.
    1388   // If this guess is wrong, check_free_list may be called.
    1389   //  int bufsize = 2*nelems0;
    1390 
    1391   memset(cdone_, -1, ncols0_);
    1392   memset(rdone_, -1, nrows0_);
    1393 
    1394   rowduals_ = new double[nrows0_];
    1395   ClpDisjointCopyN(si->getRowPrice(), nrows1, rowduals_);
    1396 
    1397   rcosts_ = new double[ncols0_];
    1398   ClpDisjointCopyN(si->getReducedCost(), ncols1, rcosts_);
    1399   if (maxmin_<0.0) {
    1400     // change so will look as if minimize
    1401     int i;
    1402     for (i=0;i<nrows1;i++)
    1403       rowduals_[i] = - rowduals_[i];
    1404     for (i=0;i<ncols1;i++) {
    1405       rcosts_[i] = - rcosts_[i];
    1406     }
    1407   }
    1408 
    1409   //ClpDisjointCopyN(si->getRowUpper(), nrows1, rup_);
    1410   //ClpDisjointCopyN(si->getRowLower(), nrows1, rlo_);
    1411 
    1412   ClpDisjointCopyN(si->getColSolution(), ncols1, sol_);
    1413   si->setDblParam(ClpObjOffset,originalOffset_);
    1414 
    1415   for (int j=0; j<ncols1; j++) {
    1416     CoinBigIndex kcs = mcstrt_[j];
    1417     CoinBigIndex kce = kcs + hincol_[j];
    1418     for (CoinBigIndex k=kcs; k<kce; ++k) {
    1419       link_[k] = k+1;
    1420     }
    1421   }
    1422   {
    1423     int ml = maxlink_;
    1424     for (CoinBigIndex k=nelemsr; k<ml; ++k)
    1425       link_[k] = k+1;
    1426     link_[ml-1] = NO_LINK;
    1427   }
    1428   free_list_ = nelemsr;
    1429 }
    1430 
    1431 
     1529
Note: See TracChangeset for help on using the changeset viewer.