source: branches/devel-1/Presolve.cpp @ 41

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

Testing presolve on nasty burlington problem

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 23.7 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 "Presolve.hpp"
17#include "PresolveMatrix.hpp"
18
19#include "PresolveEmpty.hpp"
20#include "PresolveFixed.hpp"
21#include "PresolvePsdebug.hpp"
22#include "PresolveSingleton.hpp"
23#include "PresolveDoubleton.hpp"
24#include "PresolveZeros.hpp"
25#include "PresolveSubst.hpp"
26#include "PresolveForcing.hpp"
27#include "PresolveDual.hpp"
28#include "PresolveTighten.hpp"
29#include "PresolveUseless.hpp"
30#include "PresolveDupcol.hpp"
31#include "PresolveImpliedFree.hpp"
32#include "PresolveIsolated.hpp"
33#include "ClpMessage.hpp"
34
35
36
37Presolve::Presolve() :
38  originalModel_(NULL),
39  presolvedModel_(NULL),
40  originalColumn_(NULL),
41  paction_(0),
42  ncols_(0),
43  nrows_(0),
44  nelems_(0),
45  numberPasses_(5)
46{
47}
48
49Presolve::~Presolve()
50{
51  gutsOfDestroy();
52}
53// Gets rid of presolve actions (e.g.when infeasible)
54void 
55Presolve::gutsOfDestroy()
56{
57 const PresolveAction *paction = paction_;
58  while (paction) {
59    const PresolveAction *next = paction->next;
60    delete paction;
61    paction = next;
62  }
63  delete [] originalColumn_;
64  paction_=NULL;
65  originalColumn_=NULL;
66}
67
68#if 0
69void Presolve::presolve(ClpSimplex& si)
70{
71  ncols_ = si.getNumCols();
72  nrows_ = si.getNumRows();
73  nelems_ = si.getNumElements();
74
75  double maxmin = si.getObjSense();
76
77  delete presolvedModel_;
78  presolvedModel_ = NULL;
79  originalModel_ = NULL;
80
81  PresolveMatrix prob(ncols_,
82                         maxmin,
83                         si,
84                         nrows_, nelems_,false);
85  prob.originalModel_ = originalModel_;
86
87  paction_ = presolve(&prob);
88
89  prob.update_model(si, nrows_, ncols_, nelems_);
90  int ncolsNow = si.getNumCols();
91  originalColumn_ = new int[ncolsNow];
92  memcpy(originalColumn_,prob.originalColumn_,ncolsNow*sizeof(int));
93}
94 
95
96// This is the toplevel postsolving routine.
97void Presolve::postsolve(ClpSimplex& si,
98                            unsigned char *colstat,
99                            unsigned char *rowstat)
100{
101  if (!si.isProvenOptimal()) {
102    throw CoinError("solution not optimal",
103                    "postsolve",
104                    "Presolve");
105  }
106
107  // this is the size of the original problem
108  const int ncols0  = ncols_;
109  const int nrows0  = nrows_;
110  const int nelems0 = nelems_;
111
112  double *acts = new double[nrows0];
113  double *sol  = new double[ncols0];
114
115  {
116    // put this in a block to make sure it gets deleted early
117    PostsolveMatrix prob(si,
118                            ncols0,
119                            nrows0,
120                            nelems0,
121                            si.getObjSense(),
122                            // end prepost
123
124                            sol, acts,
125                         colstat, rowstat);
126    prob.originalModel_ = originalModel_;
127
128    postsolve(prob);
129
130    // Postsolve does not use a standard matrix format, since row
131    // and col numbers change rapidly.
132    // The col rep is just a linked list.
133    int *hrow = new int[nelems0];
134    int *mcstrt = new int[ncols0+1];
135    double *colels = new double[nelems0];
136
137    {
138      const int *link = prob.link_;
139      int start = 0;
140     
141      for (int j=0; j<ncols0; j++) {
142        int k = prob.mcstrt_[j];
143        int n = prob.hincol_[j];
144
145        mcstrt[j] = start;
146        for (int i=0; i<n; i++) {
147          hrow[start] = prob.hrow_[k];
148          colels[start] = prob.colels_[k];
149          k = link[k];
150          start++;
151        }
152      }
153      mcstrt[ncols0] = start;
154    }
155   
156    si.loadProblem(ncols0, nrows0,
157                   mcstrt, hrow, colels,
158                   
159                   prob.clo_, prob.cup_,
160                   prob.cost_,
161                   prob.rlo_, prob.rup_);
162
163    delete[]hrow;
164    delete[]mcstrt;
165    delete[]colels;
166  }
167
168
169  si.setColSolution(sol);
170  delete[]acts;
171
172  delete[]sol;
173}
174#endif
175/* This version of presolve returns a pointer to a new presolved
176   model.  NULL if infeasible
177*/
178ClpSimplex * 
179Presolve::presolvedModel(ClpSimplex & si,
180                         double feasibilityTolerance,
181                         bool keepIntegers,
182                         int numberPasses)
183{
184  ncols_ = si.getNumCols();
185  nrows_ = si.getNumRows();
186  nelems_ = si.getNumElements();
187  numberPasses_ = numberPasses;
188
189  double maxmin = si.getObjSense();
190  originalModel_ = &si;
191  originalColumn_ = new int[ncols_];
192
193  // Check matrix
194  if (!si.clpMatrix()->allElementsInRange(&si,1.0e-20,1.0e20))
195    return NULL;
196
197  // result is 0 - okay, 1 infeasible, -1 go round again
198  int result = -1;
199
200  while (result==-1) {
201
202    // make new copy
203    delete presolvedModel_;
204    presolvedModel_ = new ClpSimplex(si);
205   
206
207    // drop integer information if wanted
208    if (!keepIntegers)
209      presolvedModel_->deleteIntegerInformation();
210
211   
212    PresolveMatrix prob(ncols_,
213                        maxmin,
214                        *presolvedModel_,
215                        nrows_, nelems_,true);
216    // make sure row solution correct
217    {
218      double *colels    = prob.colels_;
219      int *hrow         = prob.hrow_;
220      int *mcstrt               = prob.mcstrt_;
221      int *hincol               = prob.hincol_;
222      int ncols         = prob.ncols_;
223     
224     
225      double * csol = prob.sol_;
226      double * acts = prob.acts_;
227      int nrows = prob.nrows_;
228
229      int colx;
230
231      memset(acts,0,nrows*sizeof(double));
232     
233      for (colx = 0; colx < ncols; ++colx) {
234        double solutionValue = csol[colx];
235        for (int i=mcstrt[colx]; i<mcstrt[colx]+hincol[colx]; ++i) {
236          int row = hrow[i];
237          double coeff = colels[i];
238          acts[row] += solutionValue*coeff;
239        }
240      }
241    }
242    prob.handler_ = presolvedModel_->messageHandler();
243    prob.messages_ = presolvedModel_->messages();
244
245    // move across feasibility tolerance
246    prob.feasibilityTolerance_ = feasibilityTolerance;
247
248    // Do presolve
249
250    paction_ = presolve(&prob);
251
252    result =0; 
253
254    if (prob.status_==0&&paction_) {
255      // feasible
256   
257      prob.update_model(*presolvedModel_, nrows_, ncols_, nelems_);
258      // copy status and solution
259      memcpy(presolvedModel_->primalColumnSolution(),
260             prob.sol_,prob.ncols_*sizeof(double));
261      memcpy(presolvedModel_->primalRowSolution(),
262             prob.acts_,prob.nrows_*sizeof(double));
263      memcpy(presolvedModel_->statusArray(),
264             prob.colstat_,prob.ncols_*sizeof(unsigned char));
265      memcpy(presolvedModel_->statusArray()+prob.ncols_,
266             prob.rowstat_,prob.nrows_*sizeof(unsigned char));
267      delete [] prob.sol_;
268      delete [] prob.acts_;
269      delete [] prob.colstat_;
270      prob.sol_=NULL;
271      prob.acts_=NULL;
272      prob.colstat_=NULL;
273     
274      int ncolsNow = presolvedModel_->getNumCols();
275      memcpy(originalColumn_,prob.originalColumn_,ncolsNow*sizeof(int));
276      // now clean up integer variables.  This can modify original
277      int i;
278      const char * information = presolvedModel_->integerInformation();
279      if (information) {
280        int numberChanges=0;
281        double * lower0 = originalModel_->columnLower();
282        double * upper0 = originalModel_->columnUpper();
283        double * lower = presolvedModel_->columnLower();
284        double * upper = presolvedModel_->columnUpper();
285        for (i=0;i<ncolsNow;i++) {
286          if (!information[i])
287            continue;
288          int iOriginal = originalColumn_[i];
289          double lowerValue0 = lower0[iOriginal];
290          double upperValue0 = upper0[iOriginal];
291          double lowerValue = ceil(lower[i]-1.0e-5);
292          double upperValue = floor(upper[i]+1.0e-5);
293          lower[i]=lowerValue;
294          upper[i]=upperValue;
295          if (lowerValue>upperValue) {
296            numberChanges++;
297            presolvedModel_->messageHandler()->message(CLP_PRESOLVE_COLINFEAS,
298                                                       presolvedModel_->messages())
299                                                         <<iOriginal
300                                                         <<lowerValue
301                                                         <<upperValue
302                                                         <<CoinMessageEol;
303            result=1;
304          } else {
305            if (lowerValue>lowerValue0+1.0e-8) {
306              lower0[iOriginal] = lowerValue;
307              numberChanges++;
308            }
309            if (upperValue<upperValue0-1.0e-8) {
310              upper0[iOriginal] = upperValue;
311              numberChanges++;
312            }
313          }       
314        }
315        if (numberChanges) {
316          presolvedModel_->messageHandler()->message(CLP_PRESOLVE_INTEGERMODS,
317                                                     presolvedModel_->messages())
318                                                       <<numberChanges
319                                                       <<CoinMessageEol;
320          if (!result) {
321            result = -1; // round again
322          }
323        }
324      }
325    } else {
326      // infeasible
327      delete [] prob.sol_;
328      delete [] prob.acts_;
329      delete [] prob.colstat_;
330      result=1;
331    }
332  }
333  if (!result) {
334    int nrowsAfter = presolvedModel_->getNumRows();
335    int ncolsAfter = presolvedModel_->getNumCols();
336    int nelsAfter = presolvedModel_->getNumElements();
337    presolvedModel_->messageHandler()->message(CLP_PRESOLVE_STATS,
338                                               presolvedModel_->messages())
339                                                 <<nrowsAfter<< -(nrows_ - nrowsAfter)
340                                                 << ncolsAfter<< -(ncols_ - ncolsAfter)
341                                                 <<nelsAfter<< -(nelems_ - nelsAfter)
342                                                 <<CoinMessageEol;
343  } else {
344    gutsOfDestroy();
345    delete presolvedModel_;
346    presolvedModel_=NULL;
347  }
348  return presolvedModel_;
349}
350
351// Return pointer to presolved model
352ClpSimplex * 
353Presolve::model() const
354{
355  return presolvedModel_;
356}
357// Return pointer to original model
358ClpSimplex * 
359Presolve::originalModel() const
360{
361  return originalModel_;
362}
363void 
364Presolve::postsolve(bool updateStatus)
365{
366  if (!presolvedModel_->isProvenOptimal()) {
367    presolvedModel_->messageHandler()->message(CLP_PRESOLVE_NONOPTIMAL,
368                                             presolvedModel_->messages())
369                                               <<CoinMessageEol;
370  }
371
372  // this is the size of the original problem
373  const int ncols0  = ncols_;
374  const int nrows0  = nrows_;
375  const int nelems0 = nelems_;
376
377  // reality check
378  assert(ncols0==originalModel_->getNumCols());
379  assert(nrows0==originalModel_->getNumRows());
380
381  // this is the reduced problem
382  int ncols = presolvedModel_->getNumCols();
383  int nrows = presolvedModel_->getNumRows();
384
385  double *acts = originalModel_->primalRowSolution();
386  double *sol  = originalModel_->primalColumnSolution();
387
388  unsigned char * rowstat=NULL;
389  unsigned char * colstat = NULL;
390  if (updateStatus) {
391    unsigned char *status = originalModel_->statusArray();
392    rowstat = status + ncols0;
393    colstat = status;
394    memcpy(colstat, presolvedModel_->statusArray(), ncols);
395    memcpy(rowstat, presolvedModel_->statusArray()+ncols, nrows);
396  }
397
398
399  PostsolveMatrix prob(*presolvedModel_,
400                       ncols0,
401                       nrows0,
402                       nelems0,
403                       presolvedModel_->getObjSense(),
404                       // end prepost
405                       
406                       sol, acts,
407                       colstat, rowstat);
408   
409  postsolve(prob);
410
411}
412
413// return pointer to original columns
414const int * 
415Presolve::originalColumns() const
416{
417  return originalColumn_;
418}
419// Set pointer to original model
420void 
421Presolve::setOriginalModel(ClpSimplex * model)
422{
423  originalModel_=model;
424}
425
426// A lazy way to restrict which transformations are applied
427// during debugging.
428static int ATOI(const char *name)
429{
430  return true;
431#if     DEBUG_PRESOLVE || PRESOLVE_SUMMARY
432  if (getenv(name)) {
433    int val = atoi(getenv(name));
434    printf("%s = %d\n", name, val);
435    return (val);
436  } else {
437    if (strcmp(name,"off"))
438      return (true);
439    else
440      return (false);
441  }
442#else
443  return (true);
444#endif
445}
446//#define DEBUG_PRESOLVE 1
447#if DEBUG_PRESOLVE
448void check_sol(PresolveMatrix *prob,double tol)
449{
450  double *colels        = prob->colels_;
451  int *hrow             = prob->hrow_;
452  int *mcstrt           = prob->mcstrt_;
453  int *hincol           = prob->hincol_;
454  int *hinrow           = prob->hinrow_;
455  int ncols             = prob->ncols_;
456
457
458  double * csol = prob->sol_;
459  double * acts = prob->acts_;
460  double * clo = prob->clo_;
461  double * cup = prob->cup_;
462  int nrows = prob->nrows_;
463  double * rlo = prob->rlo_;
464  double * rup = prob->rup_;
465
466  int colx;
467
468  double * rsol = new double[nrows];
469  memset(rsol,0,nrows*sizeof(double));
470
471  for (colx = 0; colx < ncols; ++colx) {
472    if (1) {
473      int k = mcstrt[colx];
474      int nx = hincol[colx];
475      double solutionValue = csol[colx];
476      for (int i=0; i<nx; ++i) {
477        int row = hrow[k];
478        double coeff = colels[k];
479        k++;
480        rsol[row] += solutionValue*coeff;
481      }
482      if (csol[colx]<clo[colx]-tol) {
483        printf("low CSOL:  %d  - %g %g %g\n",
484                   colx, clo[colx], csol[colx], cup[colx]);
485      } else if (csol[colx]>cup[colx]+tol) {
486        printf("high CSOL:  %d  - %g %g %g\n",
487                   colx, clo[colx], csol[colx], cup[colx]);
488      } 
489    }
490  }
491  int rowx;
492  for (rowx = 0; rowx < nrows; ++rowx) {
493    if (hinrow[rowx]) {
494      if (fabs(rsol[rowx]-acts[rowx])>tol)
495        printf("inacc RSOL:  %d - %g %g (acts_ %g) %g\n",
496                   rowx,  rlo[rowx], rsol[rowx], acts[rowx], rup[rowx]);
497      if (rsol[rowx]<rlo[rowx]-tol) {
498        printf("low RSOL:  %d - %g %g %g\n",
499                   rowx,  rlo[rowx], rsol[rowx], rup[rowx]);
500      } else if (rsol[rowx]>rup[rowx]+tol ) {
501        printf("high RSOL:  %d - %g %g %g\n",
502                   rowx,  rlo[rowx], rsol[rowx], rup[rowx]);
503      } 
504    }
505  }
506  delete [] rsol;
507}
508#endif
509// This is the presolve loop.
510// It is a separate virtual function so that it can be easily
511// customized by subclassing Presolve.
512const PresolveAction *Presolve::presolve(PresolveMatrix *prob)
513{
514  paction_ = 0;
515
516  prob->status_=0; // say feasible
517
518  paction_ = make_fixed(prob, paction_);
519
520#if     CHECK_CONSISTENCY
521  presolve_links_ok(prob->rlink_, prob->mrstrt_, prob->hinrow_, prob->nrows_);
522#endif
523
524  if (!prob->status_) {
525    const bool slackd = ATOI("SLACKD");
526    //const bool forcing = ATOI("FORCING");
527    const bool doubleton = ATOI("DOUBLETON");
528    const bool forcing = ATOI("off");
529    const bool ifree = ATOI("off");
530    const bool zerocost = ATOI("off");
531    const bool dupcol = ATOI("off");
532    const bool duprow = ATOI("off");
533    const bool dual = ATOI("off");
534
535    // some things are expensive so just do once (normally)
536
537    int i;
538    // say look at all
539    for (i=0;i<nrows_;i++) 
540      prob->rowsToDo_[i]=i;
541    prob->numberRowsToDo_=nrows_;
542    for (i=0;i<ncols_;i++) 
543      prob->colsToDo_[i]=i;
544    prob->numberColsToDo_=ncols_;
545       
546
547    int iLoop;
548#if     DEBUG_PRESOLVE
549        check_sol(prob,1.0e0);
550#endif
551
552    for (iLoop=0;iLoop<numberPasses_;iLoop++) {
553#ifdef PRESOLVE_SUMMARY
554      printf("Starting major pass %d\n",iLoop+1);
555#endif
556      const PresolveAction * const paction0 = paction_;
557      // look for substitutions with no fill
558      int fill_level=2;
559      while (1) {
560        const PresolveAction * const paction1 = paction_;
561       
562        if (slackd) {
563          bool notFinished = true;
564          while (notFinished) 
565            paction_ = slack_doubleton_action::presolve(prob, paction_,
566                                                        notFinished);
567          if (prob->status_)
568            break;
569        }
570       
571        if (doubleton) {
572          paction_ = doubleton_action::presolve(prob, paction_);
573          if (prob->status_)
574            break;
575        }
576       
577        if (zerocost) {
578          paction_ = do_tighten_action::presolve(prob, paction_);
579          if (prob->status_)
580            break;
581        }
582       
583        if (forcing) {
584          paction_ = forcing_constraint_action::presolve(prob, paction_);
585          if (prob->status_)
586            break;
587        }
588       
589        if (ifree) {
590          paction_ = implied_free_action::presolve(prob, paction_,fill_level);
591          if (prob->status_)
592            break;
593        }
594       
595#if     DEBUG_PRESOLVE
596        check_sol(prob,1.0e0);
597#endif
598
599#if     CHECK_CONSISTENCY
600        presolve_links_ok(prob->rlink_, prob->mrstrt_, prob->hinrow_, 
601                          prob->nrows_);
602#endif
603       
604#if     DEBUG_PRESOLVE
605        presolve_no_zeros(prob->mcstrt_, prob->colels_, prob->hincol_, 
606                          prob->ncols_);
607#endif
608#if     CHECK_CONSISTENCY
609        prob->consistent();
610#endif
611       
612        prob->whichpass_++;
613         
614        // set up for next pass
615        // later do faster if many changes i.e. memset and memcpy
616        prob->numberRowsToDo_ = prob->numberNextRowsToDo_;
617        int kcheck;
618        bool found;
619        kcheck=-1;
620        found=false;
621        for (i=0;i<prob->numberNextRowsToDo_;i++) {
622          int index = prob->nextRowsToDo_[i];
623          prob->unsetRowChanged(index);
624          prob->rowsToDo_[i] = index;
625          if (index==kcheck) {
626            printf("row %d on list after pass %d\n",kcheck,
627                   prob->whichpass_);
628            found=true;
629          }
630        }
631        if (!found&&kcheck>=0)
632          prob->rowsToDo_[prob->numberRowsToDo_++]=kcheck;
633        prob->numberNextRowsToDo_=0;
634        prob->numberColsToDo_ = prob->numberNextColsToDo_;
635        kcheck=-1;
636        found=false;
637        for (i=0;i<prob->numberNextColsToDo_;i++) {
638          int index = prob->nextColsToDo_[i];
639          prob->unsetColChanged(index);
640          prob->colsToDo_[i] = index;
641          if (index==kcheck) {
642            printf("col %d on list after pass %d\n",kcheck,
643                   prob->whichpass_);
644            found=true;
645          }
646        }
647        if (!found&&kcheck>=0)
648          prob->colsToDo_[prob->numberColsToDo_++]=kcheck;
649        prob->numberNextColsToDo_=0;
650        if (paction_ == paction1&&fill_level>0)
651          break;
652      }
653      // say look at all
654      int i;
655      for (i=0;i<nrows_;i++) 
656        prob->rowsToDo_[i]=i;
657      prob->numberRowsToDo_=nrows_;
658      for (i=0;i<ncols_;i++) 
659        prob->colsToDo_[i]=i;
660      prob->numberColsToDo_=ncols_;
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      if (dual) {
667        int itry;
668        for (itry=0;itry<5;itry++) {
669          const PresolveAction * const paction2 = paction_;
670          paction_ = remove_dual_action::presolve(prob, paction_);
671          if (prob->status_)
672            break;
673          if (ifree) {
674            int fill_level=0; // switches off substitution
675            paction_ = implied_free_action::presolve(prob, paction_,fill_level);
676            if (prob->status_)
677              break;
678          }
679          if (paction_ == paction2)
680            break;
681        }
682      }
683
684      if (dupcol) {
685        paction_ = dupcol_action::presolve(prob, paction_);
686        if (prob->status_)
687          break;
688      }
689#if     DEBUG_PRESOLVE
690        check_sol(prob,1.0e0);
691#endif
692     
693      if (duprow) {
694        paction_ = duprow_action::presolve(prob, paction_);
695        if (prob->status_)
696          break;
697      }
698#if     DEBUG_PRESOLVE
699      check_sol(prob,1.0e0);
700#endif
701      {
702        int * hinrow = prob->hinrow_;
703        int numberDropped=0;
704        for (i=0;i<nrows_;i++) 
705          if (!hinrow[i])
706            numberDropped++;
707        printf("%d rows dropped after pass %d\n",numberDropped,
708               prob->whichpass_+1);
709      }
710      if (paction_ == paction0)
711        break;
712         
713    }
714  }
715  if (!prob->status_) {
716    paction_ = drop_zero_coefficients(prob, paction_);
717#if     DEBUG_PRESOLVE
718        check_sol(prob,1.0e0);
719#endif
720
721    paction_ = drop_empty_cols_action::presolve(prob, paction_);
722    paction_ = drop_empty_rows_action::presolve(prob, paction_);
723#if     DEBUG_PRESOLVE
724        check_sol(prob,1.0e0);
725#endif
726  }
727 
728  if (prob->status_) {
729    if (prob->status_==1)
730          presolvedModel_->messageHandler()->message(CLP_PRESOLVE_INFEAS,
731                                             presolvedModel_->messages())
732                                               <<prob->feasibilityTolerance_
733                                               <<CoinMessageEol;
734    else if (prob->status_==2)
735          presolvedModel_->messageHandler()->message(CLP_PRESOLVE_UNBOUND,
736                                             presolvedModel_->messages())
737                                               <<CoinMessageEol;
738    else
739          presolvedModel_->messageHandler()->message(CLP_PRESOLVE_INFEASUNBOUND,
740                                             presolvedModel_->messages())
741                                               <<CoinMessageEol;
742    // get rid of data
743    gutsOfDestroy();
744  }
745  return (paction_);
746}
747
748void check_djs(PostsolveMatrix *prob);
749
750
751// We could have implemented this by having each postsolve routine
752// directly call the next one, but this may make it easier to add debugging checks.
753void Presolve::postsolve(PostsolveMatrix &prob)
754{
755  const PresolveAction *paction = paction_;
756
757  if (prob.colstat_)
758    prob.check_nbasic();
759 
760#if     DEBUG_PRESOLVE
761  check_djs(&prob);
762#endif
763 
764 
765  while (paction) {
766#if     DEBUG_PRESOLVE
767    printf("POSTSOLVING %s\n", paction->name());
768#endif
769
770    paction->postsolve(&prob);
771   
772#if     DEBUG_PRESOLVE
773    if (prob.colstat_)
774      prob.check_nbasic();
775#endif
776    paction = paction->next;
777#if     DEBUG_PRESOLVE
778    check_djs(&prob);
779#endif
780  }   
781 
782#if     0 && DEBUG_PRESOLVE
783  for (i=0; i<ncols0; i++) {
784    if (!cdone[i]) {
785      printf("!cdone[%d]\n", i);
786      abort();
787    }   
788  }
789 
790  for (i=0; i<nrows0; i++) {
791    if (!rdone[i]) {
792      printf("!rdone[%d]\n", i);
793      abort();
794    }
795  }
796 
797 
798  for (i=0; i<ncols0; i++) {
799    if (sol[i] < -1e10 || sol[i] > 1e10)
800      printf("!!!%d %g\n", i, sol[i]);
801   
802  }
803 
804 
805#endif
806 
807#if     0 && DEBUG_PRESOLVE
808  // debug check:  make sure we ended up with same original matrix
809  {
810    int identical = 1;
811   
812    for (int i=0; i<ncols0; i++) {
813      PRESOLVEASSERT(hincol[i] == &prob->mcstrt0[i+1] - &prob->mcstrt0[i]);
814      int kcs0 = &prob->mcstrt0[i];
815      int kcs = mcstrt[i];
816      int n = hincol[i];
817      for (int k=0; k<n; k++) {
818        int k1 = presolve_find_row1(&prob->hrow0[kcs0+k], kcs, kcs+n, hrow);
819       
820        if (k1 == kcs+n) {
821          printf("ROW %d NOT IN COL %d\n", &prob->hrow0[kcs0+k], i);
822          abort();
823        }
824       
825        if (colels[k1] != &prob->dels0[kcs0+k])
826          printf("BAD COLEL[%d %d %d]:  %g\n",
827                 k1, i, &prob->hrow0[kcs0+k], colels[k1] - &prob->dels0[kcs0+k]);
828       
829        if (kcs0+k != k1)
830          identical=0;
831      }
832    }
833    printf("identical? %d\n", identical);
834  }
835#endif
836  // put back duals
837  memcpy(originalModel_->dualRowSolution(),prob.rowduals_,
838         nrows_*sizeof(double));
839  double maxmin = originalModel_->getObjSense();
840  if (maxmin<0.0) {
841    // swap signs
842    int i;
843    double * pi = originalModel_->dualRowSolution();
844    for (i=0;i<nrows_;i++)
845      pi[i] = -pi[i];
846  }
847  // Now check solution
848  memcpy(originalModel_->dualColumnSolution(),
849         originalModel_->objective(),ncols_*sizeof(double));
850  originalModel_->transposeTimes(-1.0,
851                                 originalModel_->dualRowSolution(),
852                                 originalModel_->dualColumnSolution());
853  memset(originalModel_->primalRowSolution(),0,nrows_*sizeof(double));
854  originalModel_->times(1.0,originalModel_->primalColumnSolution(),
855                        originalModel_->primalRowSolution());
856  originalModel_->checkSolution();
857  presolvedModel_->messageHandler()->message(CLP_PRESOLVE_POSTSOLVE,
858                                            presolvedModel_->messages())
859                                              <<originalModel_->objectiveValue()
860                                              <<originalModel_->sumDualInfeasibilities()
861                                              <<originalModel_->numberDualInfeasibilities()
862                                              <<originalModel_->sumPrimalInfeasibilities()
863                                              <<originalModel_->numberPrimalInfeasibilities()
864                                               <<CoinMessageEol;
865 
866  //originalModel_->objectiveValue_=objectiveValue_;
867  originalModel_->setNumberIterations(presolvedModel_->numberIterations());
868  if (!presolvedModel_->status()) {
869    if (!originalModel_->numberDualInfeasibilities()&&
870        !originalModel_->numberPrimalInfeasibilities()) {
871      originalModel_->setProblemStatus( 0);
872    } else {
873      originalModel_->setProblemStatus( -1);
874      presolvedModel_->messageHandler()->message(CLP_PRESOLVE_NEEDS_CLEANING,
875                                            presolvedModel_->messages())
876                                              <<CoinMessageEol;
877    }
878  } else {
879    originalModel_->setProblemStatus( presolvedModel_->status());
880  }
881}
882
883
884
885
886
887
888
889
890#if     DEBUG_PRESOLVE
891void check_djs(PostsolveMatrix *prob)
892{
893  return;
894  double *colels        = prob->colels_;
895  int *hrow             = prob->hrow_;
896  int *mcstrt           = prob->mcstrt_;
897  int *hincol           = prob->hincol_;
898  int *link             = prob->link_;
899  int ncols             = prob->ncols_;
900
901  double *dcost = prob->cost_;
902
903  double *rcosts        = prob->rcosts_;
904
905  double *rowduals = prob->rowduals_;
906
907  const double maxmin   = prob->maxmin_;
908
909  char *cdone   = prob->cdone_;
910
911  double * csol = prob->sol_;
912  double * clo = prob->clo_;
913  double * cup = prob->cup_;
914  int nrows = prob->nrows_;
915  double * rlo = prob->rlo_;
916  double * rup = prob->rup_;
917  char *rdone   = prob->rdone_;
918
919  int colx;
920
921  double * rsol = new double[nrows];
922  memset(rsol,0,nrows*sizeof(double));
923
924  for (colx = 0; colx < ncols; ++colx) {
925    if (cdone[colx]) {
926      int k = mcstrt[colx];
927      int nx = hincol[colx];
928      double dj = maxmin * dcost[colx];
929      double solutionValue = csol[colx];
930      for (int i=0; i<nx; ++i) {
931        int row = hrow[k];
932        double coeff = colels[k];
933        k = link[k];
934       
935        dj -= rowduals[row] * coeff;
936        rsol[row] += solutionValue*coeff;
937      }
938      if (! (fabs(rcosts[colx] - dj) < 100*ZTOLDP))
939        printf("BAD DJ:  %d %g %g\n",
940               colx, rcosts[colx], dj);
941      if (cup[colx]-clo[colx]>1.0e-6) {
942        if (csol[colx]<clo[colx]+1.0e-6) {
943          if (dj <-1.0e-6)
944            printf("neg DJ:  %d %g  - %g %g %g\n",
945                   colx, dj, clo[colx], csol[colx], cup[colx]);
946        } else if (csol[colx]>cup[colx]-1.0e-6) {
947          if (dj > 1.0e-6)
948            printf("pos DJ:  %d %g  - %g %g %g\n",
949                   colx, dj, clo[colx], csol[colx], cup[colx]);
950        } else {
951          if (fabs(dj) >1.0e-6)
952            printf("nonzero DJ:  %d %g  - %g %g %g\n",
953                   colx, dj, clo[colx], csol[colx], cup[colx]);
954        }
955      } 
956    }
957  }
958  int rowx;
959  for (rowx = 0; rowx < nrows; ++rowx) {
960    if (rdone[rowx]) {
961      if (rup[rowx]-rlo[rowx]>1.0e-6) {
962        double dj = rowduals[rowx];
963        if (rsol[rowx]<rlo[rowx]+1.0e-6) {
964          if (dj <-1.0e-6)
965            printf("neg rDJ:  %d %g  - %g %g %g\n",
966                   rowx, dj, rlo[rowx], rsol[rowx], rup[rowx]);
967        } else if (rsol[rowx]>rup[rowx]-1.0e-6) {
968          if (dj > 1.0e-6)
969            printf("pos rDJ:  %d %g  - %g %g %g\n",
970                   rowx, dj, rlo[rowx], rsol[rowx], rup[rowx]);
971        } else {
972          if (fabs(dj) >1.0e-6)
973            printf("nonzero rDJ:  %d %g  - %g %g %g\n",
974                   rowx, dj, rlo[rowx], rsol[rowx], rup[rowx]);
975        }
976      } 
977    }
978  }
979  delete [] rsol;
980}
981#endif
982
983
Note: See TracBrowser for help on using the repository browser.