source: trunk/ClpPresolve.cpp @ 303

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

Yet more gub

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