source: trunk/ClpPresolve.cpp @ 259

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

Hopefully improving primal and presolve

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