source: trunk/ClpPresolve.cpp @ 353

Last change on this file since 353 was 353, checked in by forrest, 16 years ago

minor

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