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

Last change on this file since 40 was 40, checked in by forrest, 18 years ago

Still trying to harden presolve

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 22.4 KB
Line 
1// Copyright (C) 2002, International Business Machines
2// Corporation and others.  All Rights Reserved.
3
4#define CHECK_CONSISTENCY       1
5
6#include <stdio.h>
7
8#include <assert.h>
9#include <iostream>
10
11#include "CoinHelperFunctions.hpp"
12
13#include "CoinPackedMatrix.hpp"
14#include "ClpSimplex.hpp"
15
16#include "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    // drop integer information if wanted
207    if (!keepIntegers)
208      presolvedModel_->deleteIntegerInformation();
209
210   
211    PresolveMatrix prob(ncols_,
212                        maxmin,
213                        *presolvedModel_,
214                        nrows_, nelems_,true);
215    prob.handler_ = presolvedModel_->messageHandler();
216    prob.messages_ = presolvedModel_->messages();
217
218    // move across feasibility tolerance
219    prob.feasibilityTolerance_ = feasibilityTolerance;
220
221    // Do presolve
222
223    paction_ = presolve(&prob);
224
225    result =0; 
226
227    if (prob.status_==0&&paction_) {
228      // feasible
229   
230      prob.update_model(*presolvedModel_, nrows_, ncols_, nelems_);
231      // copy status and solution
232      memcpy(presolvedModel_->primalColumnSolution(),
233             prob.sol_,prob.ncols_*sizeof(double));
234      memcpy(presolvedModel_->primalRowSolution(),
235             prob.acts_,prob.nrows_*sizeof(double));
236      memcpy(presolvedModel_->statusArray(),
237             prob.colstat_,prob.ncols_*sizeof(unsigned char));
238      memcpy(presolvedModel_->statusArray()+prob.ncols_,
239             prob.rowstat_,prob.nrows_*sizeof(unsigned char));
240      delete [] prob.sol_;
241      delete [] prob.acts_;
242      delete [] prob.colstat_;
243      prob.sol_=NULL;
244      prob.acts_=NULL;
245      prob.colstat_=NULL;
246     
247      int ncolsNow = presolvedModel_->getNumCols();
248      memcpy(originalColumn_,prob.originalColumn_,ncolsNow*sizeof(int));
249      // now clean up integer variables.  This can modify original
250      int i;
251      const char * information = presolvedModel_->integerInformation();
252      if (information) {
253        int numberChanges=0;
254        double * lower0 = originalModel_->columnLower();
255        double * upper0 = originalModel_->columnUpper();
256        double * lower = presolvedModel_->columnLower();
257        double * upper = presolvedModel_->columnUpper();
258        for (i=0;i<ncolsNow;i++) {
259          if (!information[i])
260            continue;
261          int iOriginal = originalColumn_[i];
262          double lowerValue0 = lower0[iOriginal];
263          double upperValue0 = upper0[iOriginal];
264          double lowerValue = ceil(lower[i]-1.0e-5);
265          double upperValue = floor(upper[i]+1.0e-5);
266          lower[i]=lowerValue;
267          upper[i]=upperValue;
268          if (lowerValue>upperValue) {
269            numberChanges++;
270            presolvedModel_->messageHandler()->message(CLP_PRESOLVE_COLINFEAS,
271                                                       presolvedModel_->messages())
272                                                         <<iOriginal
273                                                         <<lowerValue
274                                                         <<upperValue
275                                                         <<CoinMessageEol;
276            result=1;
277          } else {
278            if (lowerValue>lowerValue0+1.0e-8) {
279              lower0[iOriginal] = lowerValue;
280              numberChanges++;
281            }
282            if (upperValue<upperValue0-1.0e-8) {
283              upper0[iOriginal] = upperValue;
284              numberChanges++;
285            }
286          }       
287        }
288        if (numberChanges) {
289          presolvedModel_->messageHandler()->message(CLP_PRESOLVE_INTEGERMODS,
290                                                     presolvedModel_->messages())
291                                                       <<numberChanges
292                                                       <<CoinMessageEol;
293          if (!result) {
294            result = -1; // round again
295          }
296        }
297      }
298    } else {
299      // infeasible
300      delete [] prob.sol_;
301      delete [] prob.acts_;
302      delete [] prob.colstat_;
303      result=1;
304    }
305  }
306  if (!result) {
307    int nrowsAfter = presolvedModel_->getNumRows();
308    int ncolsAfter = presolvedModel_->getNumCols();
309    int nelsAfter = presolvedModel_->getNumElements();
310    presolvedModel_->messageHandler()->message(CLP_PRESOLVE_STATS,
311                                               presolvedModel_->messages())
312                                                 <<nrowsAfter<< -(nrows_ - nrowsAfter)
313                                                 << ncolsAfter<< -(ncols_ - ncolsAfter)
314                                                 <<nelsAfter<< -(nelems_ - nelsAfter)
315                                                 <<CoinMessageEol;
316  } else {
317    gutsOfDestroy();
318    delete presolvedModel_;
319    presolvedModel_=NULL;
320  }
321  return presolvedModel_;
322}
323
324// Return pointer to presolved model
325ClpSimplex * 
326Presolve::model() const
327{
328  return presolvedModel_;
329}
330// Return pointer to original model
331ClpSimplex * 
332Presolve::originalModel() const
333{
334  return originalModel_;
335}
336void 
337Presolve::postsolve(bool updateStatus)
338{
339  if (!presolvedModel_->isProvenOptimal()) {
340    presolvedModel_->messageHandler()->message(CLP_PRESOLVE_NONOPTIMAL,
341                                             presolvedModel_->messages())
342                                               <<CoinMessageEol;
343  }
344
345  // this is the size of the original problem
346  const int ncols0  = ncols_;
347  const int nrows0  = nrows_;
348  const int nelems0 = nelems_;
349
350  // reality check
351  assert(ncols0==originalModel_->getNumCols());
352  assert(nrows0==originalModel_->getNumRows());
353
354  // this is the reduced problem
355  int ncols = presolvedModel_->getNumCols();
356  int nrows = presolvedModel_->getNumRows();
357
358  double *acts = originalModel_->primalRowSolution();
359  double *sol  = originalModel_->primalColumnSolution();
360
361  unsigned char * rowstat=NULL;
362  unsigned char * colstat = NULL;
363  if (updateStatus) {
364    unsigned char *status = originalModel_->statusArray();
365    rowstat = status + ncols0;
366    colstat = status;
367    memcpy(colstat, presolvedModel_->statusArray(), ncols);
368    memcpy(rowstat, presolvedModel_->statusArray()+ncols, nrows);
369  }
370
371
372  PostsolveMatrix prob(*presolvedModel_,
373                       ncols0,
374                       nrows0,
375                       nelems0,
376                       presolvedModel_->getObjSense(),
377                       // end prepost
378                       
379                       sol, acts,
380                       colstat, rowstat);
381   
382  postsolve(prob);
383
384}
385
386// return pointer to original columns
387const int * 
388Presolve::originalColumns() const
389{
390  return originalColumn_;
391}
392// Set pointer to original model
393void 
394Presolve::setOriginalModel(ClpSimplex * model)
395{
396  originalModel_=model;
397}
398
399// A lazy way to restrict which transformations are applied
400// during debugging.
401static int ATOI(const char *name)
402{
403  return true;
404#if     DEBUG_PRESOLVE || PRESOLVE_SUMMARY
405  if (getenv(name)) {
406    int val = atoi(getenv(name));
407    printf("%s = %d\n", name, val);
408    return (val);
409  } else {
410    if (strcmp(name,"off"))
411      return (true);
412    else
413      return (false);
414  }
415#else
416  return (true);
417#endif
418}
419
420#ifdef  DEBUG_PRESOLVE
421void check_sol(PresolveMatrix *prob)
422{
423  double *colels        = prob->colels_;
424  int *hrow             = prob->hrow_;
425  int *mcstrt           = prob->mcstrt_;
426  int *hincol           = prob->hincol_;
427  int ncols             = prob->ncols_;
428
429
430  double * csol = prob->sol_;
431  double * clo = prob->clo_;
432  double * cup = prob->cup_;
433  int nrows = prob->nrows_;
434  double * rlo = prob->rlo_;
435  double * rup = prob->rup_;
436
437  int colx;
438
439  double * rsol = new double[nrows];
440  memset(rsol,0,nrows*sizeof(double));
441
442  for (colx = 0; colx < ncols; ++colx) {
443    if (1) {
444      int k = mcstrt[colx];
445      int nx = hincol[colx];
446      double solutionValue = csol[colx];
447      for (int i=0; i<nx; ++i) {
448        int row = hrow[k];
449        double coeff = colels[k];
450        k++;
451        rsol[row] += solutionValue*coeff;
452      }
453      if (csol[colx]<clo[colx]-1.0e-4) {
454        printf("low CSOL:  %d  - %g %g %g\n",
455                   colx, clo[colx], csol[colx], cup[colx]);
456      } else if (csol[colx]>cup[colx]+1.0e-4) {
457        printf("high CSOL:  %d  - %g %g %g\n",
458                   colx, clo[colx], csol[colx], cup[colx]);
459      } 
460    }
461  }
462  int rowx;
463  for (rowx = 0; rowx < nrows; ++rowx) {
464    if (1) {
465      if (rsol[rowx]<rlo[rowx]-1.0e-4) {
466        printf("low RSOL:  %d - %g %g %g\n",
467                   rowx,  rlo[rowx], rsol[rowx], rup[rowx]);
468      } else if (rsol[rowx]>rup[rowx]+1.0e-4) {
469        printf("high RSOL:  %d - %g %g %g\n",
470                   rowx,  rlo[rowx], rsol[rowx], rup[rowx]);
471      } 
472    }
473  }
474  delete [] rsol;
475}
476#endif
477// This is the presolve loop.
478// It is a separate virtual function so that it can be easily
479// customized by subclassing Presolve.
480const PresolveAction *Presolve::presolve(PresolveMatrix *prob)
481{
482  paction_ = 0;
483
484  prob->status_=0; // say feasible
485
486  paction_ = make_fixed(prob, paction_);
487
488#if     CHECK_CONSISTENCY
489  presolve_links_ok(prob->rlink_, prob->mrstrt_, prob->hinrow_, prob->nrows_);
490#endif
491
492  if (!prob->status_) {
493    const bool slackd = ATOI("SLACKD");
494    //const bool forcing = ATOI("FORCING");
495    const bool doubleton = ATOI("DOUBLETON");
496    const bool forcing = ATOI("off");
497    const bool ifree = ATOI("off");
498    const bool zerocost = ATOI("off");
499    const bool dupcol = ATOI("off");
500    const bool duprow = ATOI("off");
501    const bool dual = ATOI("off");
502
503    // some things are expensive so just do once (normally)
504
505    int i;
506    // say look at all
507    for (i=0;i<nrows_;i++) 
508      prob->rowsToDo_[i]=i;
509    prob->numberRowsToDo_=nrows_;
510    for (i=0;i<ncols_;i++) 
511      prob->colsToDo_[i]=i;
512    prob->numberColsToDo_=ncols_;
513       
514
515    int iLoop;
516
517    for (iLoop=0;iLoop<numberPasses_;iLoop++) {
518#ifdef PRESOLVE_SUMMARY
519      printf("Starting major pass %d\n",iLoop+1);
520#endif
521      const PresolveAction * const paction0 = paction_;
522      // look for substitutions with no fill
523      int fill_level=2;
524      while (1) {
525        const PresolveAction * const paction1 = paction_;
526       
527        if (slackd) {
528          bool notFinished = true;
529          while (notFinished) 
530            paction_ = slack_doubleton_action::presolve(prob, paction_,
531                                                        notFinished);
532          if (prob->status_)
533            break;
534        }
535       
536        if (doubleton) {
537          paction_ = doubleton_action::presolve(prob, paction_);
538          if (prob->status_)
539            break;
540        }
541       
542        if (zerocost) {
543          paction_ = do_tighten_action::presolve(prob, paction_);
544          if (prob->status_)
545            break;
546        }
547       
548        if (forcing) {
549          paction_ = forcing_constraint_action::presolve(prob, paction_);
550          if (prob->status_)
551            break;
552        }
553       
554        if (ifree) {
555          paction_ = implied_free_action::presolve(prob, paction_,fill_level);
556          if (prob->status_)
557            break;
558        }
559       
560       
561#if     DEBUG_PRESOLVE
562        check_sol(prob);
563#endif
564
565#if     CHECK_CONSISTENCY
566        presolve_links_ok(prob->rlink_, prob->mrstrt_, prob->hinrow_, 
567                          prob->nrows_);
568#endif
569       
570#if     DEBUG_PRESOLVE
571        presolve_no_zeros(prob->mcstrt_, prob->colels_, prob->hincol_, 
572                          prob->ncols_);
573#endif
574#if     CHECK_CONSISTENCY
575        prob->consistent();
576#endif
577       
578#if     PRESOLVE_SUMMARY
579        prob->whichpass_++;
580#endif
581         
582        // set up for next pass
583        // later do faster if many changes i.e. memset and memcpy
584        prob->numberRowsToDo_ = prob->numberNextRowsToDo_;
585        int kcheck;
586        bool found;
587        kcheck=-1;
588        found=false;
589        for (i=0;i<prob->numberNextRowsToDo_;i++) {
590          int index = prob->nextRowsToDo_[i];
591          prob->unsetRowChanged(index);
592          prob->rowsToDo_[i] = index;
593          if (index==kcheck) {
594            printf("row %d on list after pass %d\n",kcheck,
595                   prob->whichpass_);
596            found=true;
597          }
598        }
599        if (!found&&kcheck>=0)
600          prob->rowsToDo_[prob->numberRowsToDo_++]=kcheck;
601        prob->numberNextRowsToDo_=0;
602        prob->numberColsToDo_ = prob->numberNextColsToDo_;
603        kcheck=-1;
604        found=false;
605        for (i=0;i<prob->numberNextColsToDo_;i++) {
606          int index = prob->nextColsToDo_[i];
607          prob->unsetColChanged(index);
608          prob->colsToDo_[i] = index;
609          if (index==kcheck) {
610            printf("col %d on list after pass %d\n",kcheck,
611                   prob->whichpass_);
612            found=true;
613          }
614        }
615        if (!found&&kcheck>=0)
616          prob->colsToDo_[prob->numberColsToDo_++]=kcheck;
617        prob->numberNextColsToDo_=0;
618        if (paction_ == paction1&&fill_level>0)
619          break;
620      }
621      // say look at all
622      int i;
623      for (i=0;i<nrows_;i++) 
624        prob->rowsToDo_[i]=i;
625      prob->numberRowsToDo_=nrows_;
626      for (i=0;i<ncols_;i++) 
627        prob->colsToDo_[i]=i;
628      prob->numberColsToDo_=ncols_;
629      // now expensive things
630      // this caused world.mps to run into numerical difficulties
631#ifdef PRESOLVE_SUMMARY
632      printf("Starting expensive\n");
633#endif
634
635      if (dual) {
636        int itry;
637        for (itry=0;itry<5;itry++) {
638          const PresolveAction * const paction2 = paction_;
639          paction_ = remove_dual_action::presolve(prob, paction_);
640          if (prob->status_)
641            break;
642          if (ifree) {
643            int fill_level=0; // switches off substitution
644            paction_ = implied_free_action::presolve(prob, paction_,fill_level);
645            if (prob->status_)
646              break;
647          }
648          if (paction_ == paction2)
649            break;
650        }
651      }
652
653      if (dupcol) {
654        paction_ = dupcol_action::presolve(prob, paction_);
655        if (prob->status_)
656          break;
657      }
658     
659      if (duprow) {
660        paction_ = duprow_action::presolve(prob, paction_);
661        if (prob->status_)
662          break;
663      }
664      if (paction_ == paction0)
665        break;
666         
667    }
668  }
669  if (!prob->status_) {
670    paction_ = drop_zero_coefficients(prob, paction_);
671
672    paction_ = drop_empty_cols_action::presolve(prob, paction_);
673    paction_ = drop_empty_rows_action::presolve(prob, paction_);
674  }
675 
676  if (prob->status_) {
677    if (prob->status_==1)
678          presolvedModel_->messageHandler()->message(CLP_PRESOLVE_INFEAS,
679                                             presolvedModel_->messages())
680                                               <<prob->feasibilityTolerance_
681                                               <<CoinMessageEol;
682    else if (prob->status_==2)
683          presolvedModel_->messageHandler()->message(CLP_PRESOLVE_UNBOUND,
684                                             presolvedModel_->messages())
685                                               <<CoinMessageEol;
686    else
687          presolvedModel_->messageHandler()->message(CLP_PRESOLVE_INFEASUNBOUND,
688                                             presolvedModel_->messages())
689                                               <<CoinMessageEol;
690    // get rid of data
691    gutsOfDestroy();
692  }
693  return (paction_);
694}
695
696void check_djs(PostsolveMatrix *prob);
697
698
699// We could have implemented this by having each postsolve routine
700// directly call the next one, but this may make it easier to add debugging checks.
701void Presolve::postsolve(PostsolveMatrix &prob)
702{
703  const PresolveAction *paction = paction_;
704
705  if (prob.colstat_)
706    prob.check_nbasic();
707 
708#if     DEBUG_PRESOLVE
709  check_djs(&prob);
710#endif
711 
712 
713  while (paction) {
714#if     DEBUG_PRESOLVE
715    printf("POSTSOLVING %s\n", paction->name());
716#endif
717
718    paction->postsolve(&prob);
719   
720#if     DEBUG_PRESOLVE
721    if (prob.colstat_)
722      prob.check_nbasic();
723#endif
724    paction = paction->next;
725#if     DEBUG_PRESOLVE
726    check_djs(&prob);
727#endif
728  }   
729 
730#if     0 && DEBUG_PRESOLVE
731  for (i=0; i<ncols0; i++) {
732    if (!cdone[i]) {
733      printf("!cdone[%d]\n", i);
734      abort();
735    }   
736  }
737 
738  for (i=0; i<nrows0; i++) {
739    if (!rdone[i]) {
740      printf("!rdone[%d]\n", i);
741      abort();
742    }
743  }
744 
745 
746  for (i=0; i<ncols0; i++) {
747    if (sol[i] < -1e10 || sol[i] > 1e10)
748      printf("!!!%d %g\n", i, sol[i]);
749   
750  }
751 
752 
753#endif
754 
755#if     0 && DEBUG_PRESOLVE
756  // debug check:  make sure we ended up with same original matrix
757  {
758    int identical = 1;
759   
760    for (int i=0; i<ncols0; i++) {
761      PRESOLVEASSERT(hincol[i] == &prob->mcstrt0[i+1] - &prob->mcstrt0[i]);
762      int kcs0 = &prob->mcstrt0[i];
763      int kcs = mcstrt[i];
764      int n = hincol[i];
765      for (int k=0; k<n; k++) {
766        int k1 = presolve_find_row1(&prob->hrow0[kcs0+k], kcs, kcs+n, hrow);
767       
768        if (k1 == kcs+n) {
769          printf("ROW %d NOT IN COL %d\n", &prob->hrow0[kcs0+k], i);
770          abort();
771        }
772       
773        if (colels[k1] != &prob->dels0[kcs0+k])
774          printf("BAD COLEL[%d %d %d]:  %g\n",
775                 k1, i, &prob->hrow0[kcs0+k], colels[k1] - &prob->dels0[kcs0+k]);
776       
777        if (kcs0+k != k1)
778          identical=0;
779      }
780    }
781    printf("identical? %d\n", identical);
782  }
783#endif
784  // put back duals
785  memcpy(originalModel_->dualRowSolution(),prob.rowduals_,
786         nrows_*sizeof(double));
787  double maxmin = originalModel_->getObjSense();
788  if (maxmin<0.0) {
789    // swap signs
790    int i;
791    double * pi = originalModel_->dualRowSolution();
792    for (i=0;i<nrows_;i++)
793      pi[i] = -pi[i];
794  }
795  // Now check solution
796  memcpy(originalModel_->dualColumnSolution(),
797         originalModel_->objective(),ncols_*sizeof(double));
798  originalModel_->transposeTimes(-1.0,
799                                 originalModel_->dualRowSolution(),
800                                 originalModel_->dualColumnSolution());
801  memset(originalModel_->primalRowSolution(),0,nrows_*sizeof(double));
802  originalModel_->times(1.0,originalModel_->primalColumnSolution(),
803                        originalModel_->primalRowSolution());
804  originalModel_->checkSolution();
805  presolvedModel_->messageHandler()->message(CLP_PRESOLVE_POSTSOLVE,
806                                            presolvedModel_->messages())
807                                              <<originalModel_->objectiveValue()
808                                              <<originalModel_->sumDualInfeasibilities()
809                                              <<originalModel_->numberDualInfeasibilities()
810                                              <<originalModel_->sumPrimalInfeasibilities()
811                                              <<originalModel_->numberPrimalInfeasibilities()
812                                               <<CoinMessageEol;
813 
814  //originalModel_->objectiveValue_=objectiveValue_;
815  originalModel_->setNumberIterations(presolvedModel_->numberIterations());
816  if (!presolvedModel_->status()) {
817    if (!originalModel_->numberDualInfeasibilities()&&
818        !originalModel_->numberPrimalInfeasibilities()) {
819      originalModel_->setProblemStatus( 0);
820    } else {
821      originalModel_->setProblemStatus( -1);
822      presolvedModel_->messageHandler()->message(CLP_PRESOLVE_NEEDS_CLEANING,
823                                            presolvedModel_->messages())
824                                              <<CoinMessageEol;
825    }
826  } else {
827    originalModel_->setProblemStatus( presolvedModel_->status());
828  }
829}
830
831
832
833
834
835
836
837
838#if     DEBUG_PRESOLVE
839void check_djs(PostsolveMatrix *prob)
840{
841  double *colels        = prob->colels_;
842  int *hrow             = prob->hrow_;
843  int *mcstrt           = prob->mcstrt_;
844  int *hincol           = prob->hincol_;
845  int *link             = prob->link_;
846  int ncols             = prob->ncols_;
847
848  double *dcost = prob->cost_;
849
850  double *rcosts        = prob->rcosts_;
851
852  double *rowduals = prob->rowduals_;
853
854  const double maxmin   = prob->maxmin_;
855
856  char *cdone   = prob->cdone_;
857
858  double * csol = prob->sol_;
859  double * clo = prob->clo_;
860  double * cup = prob->cup_;
861  int nrows = prob->nrows_;
862  double * rlo = prob->rlo_;
863  double * rup = prob->rup_;
864  char *rdone   = prob->rdone_;
865
866  int colx;
867
868  double * rsol = new double[nrows];
869  memset(rsol,0,nrows*sizeof(double));
870
871  for (colx = 0; colx < ncols; ++colx) {
872    if (cdone[colx]) {
873      int k = mcstrt[colx];
874      int nx = hincol[colx];
875      double dj = maxmin * dcost[colx];
876      double solutionValue = csol[colx];
877      for (int i=0; i<nx; ++i) {
878        int row = hrow[k];
879        double coeff = colels[k];
880        k = link[k];
881       
882        dj -= rowduals[row] * coeff;
883        rsol[row] += solutionValue*coeff;
884      }
885      if (! (fabs(rcosts[colx] - dj) < 100*ZTOLDP))
886        printf("BAD DJ:  %d %g %g\n",
887               colx, rcosts[colx], dj);
888      if (cup[colx]-clo[colx]>1.0e-6) {
889        if (csol[colx]<clo[colx]+1.0e-6) {
890          if (dj <-1.0e-6)
891            printf("neg DJ:  %d %g  - %g %g %g\n",
892                   colx, dj, clo[colx], csol[colx], cup[colx]);
893        } else if (csol[colx]>cup[colx]-1.0e-6) {
894          if (dj > 1.0e-6)
895            printf("pos DJ:  %d %g  - %g %g %g\n",
896                   colx, dj, clo[colx], csol[colx], cup[colx]);
897        } else {
898          if (fabs(dj) >1.0e-6)
899            printf("nonzero DJ:  %d %g  - %g %g %g\n",
900                   colx, dj, clo[colx], csol[colx], cup[colx]);
901        }
902      } 
903    }
904  }
905  int rowx;
906  for (rowx = 0; rowx < nrows; ++rowx) {
907    if (rdone[rowx]) {
908      if (rup[rowx]-rlo[rowx]>1.0e-6) {
909        double dj = rowduals[rowx];
910        if (rsol[rowx]<rlo[rowx]+1.0e-6) {
911          if (dj <-1.0e-6)
912            printf("neg rDJ:  %d %g  - %g %g %g\n",
913                   rowx, dj, rlo[rowx], rsol[rowx], rup[rowx]);
914        } else if (rsol[rowx]>rup[rowx]-1.0e-6) {
915          if (dj > 1.0e-6)
916            printf("pos rDJ:  %d %g  - %g %g %g\n",
917                   rowx, dj, rlo[rowx], rsol[rowx], rup[rowx]);
918        } else {
919          if (fabs(dj) >1.0e-6)
920            printf("nonzero rDJ:  %d %g  - %g %g %g\n",
921                   rowx, dj, rlo[rowx], rsol[rowx], rup[rowx]);
922        }
923      } 
924    }
925  }
926  delete [] rsol;
927}
928#endif
929
930
Note: See TracBrowser for help on using the repository browser.