source: branches/pre/ClpPresolve.cpp @ 219

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

pdco seems to work

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