source: branches/pre/ClpPresolve.cpp @ 212

Last change on this file since 212 was 212, checked in by forrest, 17 years ago

lots of stuff

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