source: trunk/ClpPresolve.cpp @ 238

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

Messages and stuff

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 39.0 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 returns a pointer to a new presolved
92   model and save original data to file.  Returns non-zero if infeasible
93*/
94int
95ClpPresolve::presolvedModel(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    for (iLoop=0;iLoop<numberPasses_;iLoop++) {
442#ifdef PRESOLVE_SUMMARY
443      printf("Starting major pass %d\n",iLoop+1);
444#endif
445      const CoinPresolveAction * const paction0 = paction_;
446      // look for substitutions with no fill
447      int fill_level=2;
448      //fill_level=10;
449      //printf("** fill_level == 10 !\n");
450      int whichPass=0;
451      while (1) {
452        whichPass++;
453        const CoinPresolveAction * const paction1 = paction_;
454
455        if (slackd) {
456          bool notFinished = true;
457          while (notFinished) 
458            paction_ = slack_doubleton_action::presolve(prob, paction_,
459                                                        notFinished);
460          if (prob->status_)
461            break;
462        }
463
464        if (doubleton) {
465          paction_ = doubleton_action::presolve(prob, paction_);
466          if (prob->status_)
467            break;
468        }
469
470        if (tripleton) {
471          paction_ = tripleton_action::presolve(prob, paction_);
472          if (prob->status_)
473            break;
474        }
475
476        if (zerocost) {
477          paction_ = do_tighten_action::presolve(prob, paction_);
478          if (prob->status_)
479            break;
480        }
481
482        if (forcing) {
483          paction_ = forcing_constraint_action::presolve(prob, paction_);
484          if (prob->status_)
485            break;
486        }
487
488        if (ifree) {
489          paction_ = implied_free_action::presolve(prob, paction_,fill_level);
490          if (prob->status_)
491            break;
492        }
493
494#if     DEBUG_PRESOLVE
495        check_sol(prob,1.0e0);
496#endif
497
498#if     CHECK_CONSISTENCY
499        presolve_links_ok(prob->rlink_, prob->mrstrt_, prob->hinrow_, 
500                          prob->nrows_);
501#endif
502
503#if     DEBUG_PRESOLVE
504        presolve_no_zeros(prob->mcstrt_, prob->colels_, prob->hincol_, 
505                          prob->ncols_);
506#endif
507#if     CHECK_CONSISTENCY
508        prob->consistent();
509#endif
510
511         
512        // set up for next pass
513        // later do faster if many changes i.e. memset and memcpy
514        prob->numberRowsToDo_ = prob->numberNextRowsToDo_;
515        int kcheck;
516        bool found=false;
517        kcheck=-1;
518        for (i=0;i<prob->numberNextRowsToDo_;i++) {
519          int index = prob->nextRowsToDo_[i];
520          prob->unsetRowChanged(index);
521          prob->rowsToDo_[i] = index;
522          if (index==kcheck) {
523            printf("row %d on list after pass %d\n",kcheck,
524                   whichPass);
525            found=true;
526          }
527        }
528        if (!found&&kcheck>=0)
529          prob->rowsToDo_[prob->numberRowsToDo_++]=kcheck;
530        prob->numberNextRowsToDo_=0;
531        prob->numberColsToDo_ = prob->numberNextColsToDo_;
532        kcheck=-1;
533        found=false;
534        for (i=0;i<prob->numberNextColsToDo_;i++) {
535          int index = prob->nextColsToDo_[i];
536          prob->unsetColChanged(index);
537          prob->colsToDo_[i] = index;
538          if (index==kcheck) {
539            printf("col %d on list after pass %d\n",kcheck,
540                   whichPass);
541            found=true;
542          }
543        }
544        if (!found&&kcheck>=0)
545          prob->colsToDo_[prob->numberColsToDo_++]=kcheck;
546        prob->numberNextColsToDo_=0;
547        if (paction_ == paction1&&fill_level>0)
548          break;
549      }
550      // say look at all
551      int i;
552      if (!prob->anyProhibited()) {
553        for (i=0;i<nrows_;i++) 
554          prob->rowsToDo_[i]=i;
555        prob->numberRowsToDo_=nrows_;
556        for (i=0;i<ncols_;i++) 
557          prob->colsToDo_[i]=i;
558        prob->numberColsToDo_=ncols_;
559      } else {
560        // some stuff must be left alone
561        prob->numberRowsToDo_=0;
562        for (i=0;i<nrows_;i++) 
563          if (!prob->rowProhibited(i))
564            prob->rowsToDo_[prob->numberRowsToDo_++]=i;
565        prob->numberColsToDo_=0;
566        for (i=0;i<ncols_;i++) 
567          if (!prob->colProhibited(i))
568            prob->colsToDo_[prob->numberColsToDo_++]=i;
569      }
570      // now expensive things
571      // this caused world.mps to run into numerical difficulties
572#ifdef PRESOLVE_SUMMARY
573      printf("Starting expensive\n");
574#endif
575
576      if (dual) {
577        int itry;
578        for (itry=0;itry<5;itry++) {
579          const CoinPresolveAction * const paction2 = paction_;
580          paction_ = remove_dual_action::presolve(prob, paction_);
581          if (prob->status_)
582            break;
583          if (ifree) {
584            int fill_level=0; // switches off substitution
585            paction_ = implied_free_action::presolve(prob, paction_,fill_level);
586            if (prob->status_)
587              break;
588          }
589          if (paction_ == paction2)
590            break;
591        }
592      }
593#if     DEBUG_PRESOLVE
594      check_sol(prob,1.0e0);
595#endif
596      if (dupcol) {
597        paction_ = dupcol_action::presolve(prob, paction_);
598        if (prob->status_)
599          break;
600      }
601#if     DEBUG_PRESOLVE
602        check_sol(prob,1.0e0);
603#endif
604     
605      if (duprow) {
606        paction_ = duprow_action::presolve(prob, paction_);
607        if (prob->status_)
608          break;
609      }
610#if     DEBUG_PRESOLVE
611      check_sol(prob,1.0e0);
612#endif
613      {
614        int * hinrow = prob->hinrow_;
615        int numberDropped=0;
616        for (i=0;i<nrows_;i++) 
617          if (!hinrow[i])
618            numberDropped++;
619        //printf("%d rows dropped after pass %d\n",numberDropped,
620        //     iLoop+1);
621        if (numberDropped==lastDropped)
622          break;
623        else
624          lastDropped = numberDropped;
625      }
626      if (paction_ == paction0)
627        break;
628         
629    }
630  }
631  if (!prob->status_) {
632    paction_ = drop_zero_coefficients(prob, paction_);
633#if     DEBUG_PRESOLVE
634        check_sol(prob,1.0e0);
635#endif
636
637    paction_ = drop_empty_cols_action::presolve(prob, paction_);
638    paction_ = drop_empty_rows_action::presolve(prob, paction_);
639#if     DEBUG_PRESOLVE
640        check_sol(prob,1.0e0);
641#endif
642  }
643 
644  // Messages
645  CoinMessages messages = CoinMessage(prob->messages().language());
646  if (prob->status_) {
647    if (prob->status_==1)
648          prob->messageHandler()->message(COIN_PRESOLVE_INFEAS,
649                                             messages)
650                                               <<prob->feasibilityTolerance_
651                                               <<CoinMessageEol;
652    else if (prob->status_==2)
653          prob->messageHandler()->message(COIN_PRESOLVE_UNBOUND,
654                                             messages)
655                                               <<CoinMessageEol;
656    else
657          prob->messageHandler()->message(COIN_PRESOLVE_INFEASUNBOUND,
658                                             messages)
659                                               <<CoinMessageEol;
660    // get rid of data
661    gutsOfDestroy();
662  }
663  return (paction_);
664}
665
666void check_djs(CoinPostsolveMatrix *prob);
667
668
669// We could have implemented this by having each postsolve routine
670// directly call the next one, but this may make it easier to add debugging checks.
671void ClpPresolve::postsolve(CoinPostsolveMatrix &prob)
672{
673  {
674    // Check activities
675    double *colels      = prob.colels_;
676    int *hrow           = prob.hrow_;
677    int *mcstrt         = prob.mcstrt_;
678    int *hincol         = prob.hincol_;
679    int *link           = prob.link_;
680    int ncols           = prob.ncols_;
681
682    char *cdone = prob.cdone_;
683
684    double * csol = prob.sol_;
685    int nrows = prob.nrows_;
686
687    int colx;
688
689    double * rsol = prob.acts_;
690    memset(rsol,0,nrows*sizeof(double));
691
692    for (colx = 0; colx < ncols; ++colx) {
693      if (cdone[colx]) {
694        CoinBigIndex k = mcstrt[colx];
695        int nx = hincol[colx];
696        double solutionValue = csol[colx];
697        for (int i=0; i<nx; ++i) {
698          int row = hrow[k];
699          double coeff = colels[k];
700          k = link[k];
701          rsol[row] += solutionValue*coeff;
702        }
703      }
704    }
705  }
706  const CoinPresolveAction *paction = paction_;
707
708  if (prob.colstat_)
709    prob.check_nbasic();
710 
711#if     DEBUG_PRESOLVE
712  check_djs(&prob);
713#endif
714 
715 
716  while (paction) {
717#if     DEBUG_PRESOLVE
718    printf("POSTSOLVING %s\n", paction->name());
719#endif
720
721    paction->postsolve(&prob);
722   
723#if     DEBUG_PRESOLVE
724    if (prob.colstat_)
725      prob.check_nbasic();
726#endif
727    paction = paction->next;
728#if     DEBUG_PRESOLVE
729    check_djs(&prob);
730#endif
731  }   
732 
733#if     0 && DEBUG_PRESOLVE
734  for (i=0; i<ncols0; i++) {
735    if (!cdone[i]) {
736      printf("!cdone[%d]\n", i);
737      abort();
738    }
739  }
740 
741  for (i=0; i<nrows0; i++) {
742    if (!rdone[i]) {
743      printf("!rdone[%d]\n", i);
744      abort();
745    }
746  }
747 
748 
749  for (i=0; i<ncols0; i++) {
750    if (sol[i] < -1e10 || sol[i] > 1e10)
751      printf("!!!%d %g\n", i, sol[i]);
752   
753  }
754 
755 
756#endif
757 
758#if     0 && DEBUG_PRESOLVE
759  // debug check:  make sure we ended up with same original matrix
760  {
761    int identical = 1;
762   
763    for (int i=0; i<ncols0; i++) {
764      PRESOLVEASSERT(hincol[i] == &prob->mcstrt0[i+1] - &prob->mcstrt0[i]);
765      CoinBigIndex kcs0 = &prob->mcstrt0[i];
766      CoinBigIndex kcs = mcstrt[i];
767      int n = hincol[i];
768      for (int k=0; k<n; k++) {
769        CoinBigIndex k1 = presolve_find_row1(&prob->hrow0[kcs0+k], kcs, kcs+n, hrow);
770
771        if (k1 == kcs+n) {
772          printf("ROW %d NOT IN COL %d\n", &prob->hrow0[kcs0+k], i);
773          abort();
774        }
775
776        if (colels[k1] != &prob->dels0[kcs0+k])
777          printf("BAD COLEL[%d %d %d]:  %g\n",
778                 k1, i, &prob->hrow0[kcs0+k], colels[k1] - &prob->dels0[kcs0+k]);
779
780        if (kcs0+k != k1)
781          identical=0;
782      }
783    }
784    printf("identical? %d\n", identical);
785  }
786#endif
787}
788
789
790
791
792
793
794
795
796#if     DEBUG_PRESOLVE
797void check_djs(CoinPostsolveMatrix *prob)
798{
799  //return;
800  double *colels        = prob->colels_;
801  int *hrow             = prob->hrow_;
802  int *mcstrt           = prob->mcstrt_;
803  int *hincol           = prob->hincol_;
804  int *link             = prob->link_;
805  int ncols             = prob->ncols_;
806
807  double *dcost = prob->cost_;
808
809  double *rcosts        = prob->rcosts_;
810
811  double *rowduals = prob->rowduals_;
812
813  const double maxmin   = prob->maxmin_;
814
815  char *cdone   = prob->cdone_;
816
817  double * csol = prob->sol_;
818  double * clo = prob->clo_;
819  double * cup = prob->cup_;
820  int nrows = prob->nrows_;
821  double * rlo = prob->rlo_;
822  double * rup = prob->rup_;
823  char *rdone   = prob->rdone_;
824
825  int colx;
826
827  double * rsol = new double[nrows];
828  memset(rsol,0,nrows*sizeof(double));
829
830  for (colx = 0; colx < ncols; ++colx) {
831    if (cdone[colx]) {
832      CoinBigIndex k = mcstrt[colx];
833      int nx = hincol[colx];
834      double dj = maxmin * dcost[colx];
835      double solutionValue = csol[colx];
836      for (int i=0; i<nx; ++i) {
837        int row = hrow[k];
838        double coeff = colels[k];
839        k = link[k];
840        dj -= rowduals[row] * coeff;
841        rsol[row] += solutionValue*coeff;
842      }
843      if (! (fabs(rcosts[colx] - dj) < 1.0e-4))
844        printf("BAD DJ:  %d %g %g\n",
845               colx, rcosts[colx], dj);
846      if (cup[colx]-clo[colx]>1.0e-6) {
847        if (csol[colx]<clo[colx]+1.0e-6) {
848          if (dj <-1.0e-6)
849            printf("neg DJ:  %d %g  - %g %g %g\n",
850                   colx, dj, clo[colx], csol[colx], cup[colx]);
851        } else if (csol[colx]>cup[colx]-1.0e-6) {
852          if (dj > 1.0e-6)
853            printf("pos DJ:  %d %g  - %g %g %g\n",
854                   colx, dj, clo[colx], csol[colx], cup[colx]);
855        } else {
856          if (fabs(dj) >1.0e-6)
857            printf("nonzero DJ:  %d %g  - %g %g %g\n",
858                   colx, dj, clo[colx], csol[colx], cup[colx]);
859        }
860      } 
861    }
862  }
863  int rowx;
864  for (rowx = 0; rowx < nrows; ++rowx) {
865    if (rdone[rowx]) {
866      if (rup[rowx]-rlo[rowx]>1.0e-6) {
867        double dj = rowduals[rowx];
868        if (rsol[rowx]<rlo[rowx]+1.0e-6) {
869          if (dj <-1.0e-5)
870            printf("neg rDJ:  %d %g  - %g %g %g\n",
871                   rowx, dj, rlo[rowx], rsol[rowx], rup[rowx]);
872        } else if (rsol[rowx]>rup[rowx]-1.0e-6) {
873          if (dj > 1.0e-5)
874            printf("pos rDJ:  %d %g  - %g %g %g\n",
875                   rowx, dj, rlo[rowx], rsol[rowx], rup[rowx]);
876        } else {
877          if (fabs(dj) >1.0e-5)
878            printf("nonzero rDJ:  %d %g  - %g %g %g\n",
879                   rowx, dj, rlo[rowx], rsol[rowx], rup[rowx]);
880        }
881      } 
882    }
883  }
884  delete [] rsol;
885}
886#endif
887
888static inline double getTolerance(const ClpSimplex  *si, ClpDblParam key)
889{
890  double tol;
891  if (! si->getDblParam(key, tol)) {
892    CoinPresolveAction::throwCoinError("getDblParam failed",
893                                      "CoinPrePostsolveMatrix::CoinPrePostsolveMatrix");
894  }
895  return (tol);
896}
897
898
899// Assumptions:
900// 1. nrows>=m.getNumRows()
901// 2. ncols>=m.getNumCols()
902//
903// In presolve, these values are equal.
904// In postsolve, they may be inequal, since the reduced problem
905// may be smaller, but we need room for the large problem.
906// ncols may be larger than si.getNumCols() in postsolve,
907// this at that point si will be the reduced problem,
908// but we need to reserve enough space for the original problem.
909CoinPrePostsolveMatrix::CoinPrePostsolveMatrix(const ClpSimplex * si,
910                                             int ncols_in,
911                                             int nrows_in,
912                                             CoinBigIndex nelems_in) :
913  ncols_(si->getNumCols()),
914  ncols0_(ncols_in),
915  nelems_(si->getNumElements()),
916
917  mcstrt_(new CoinBigIndex[ncols_in+1]),
918  hincol_(new int[ncols_in+1]),
919  hrow_  (new int   [2*nelems_in]),
920  colels_(new double[2*nelems_in]),
921
922  cost_(new double[ncols_in]),
923  clo_(new double[ncols_in]),
924  cup_(new double[ncols_in]),
925  rlo_(new double[nrows_in]),
926  rup_(new double[nrows_in]),
927  originalColumn_(new int[ncols_in]),
928  originalRow_(new int[nrows_in]),
929
930  ztolzb_(getTolerance(si, ClpPrimalTolerance)),
931  ztoldj_(getTolerance(si, ClpDualTolerance)),
932
933  maxmin_(si->getObjSense())
934
935{
936  si->getDblParam(ClpObjOffset,originalOffset_);
937  int ncols = si->getNumCols();
938  int nrows = si->getNumRows();
939
940  ClpDisjointCopyN(si->getColLower(), ncols, clo_);
941  ClpDisjointCopyN(si->getColUpper(), ncols, cup_);
942  ClpDisjointCopyN(si->getObjCoefficients(), ncols, cost_);
943  ClpDisjointCopyN(si->getRowLower(), nrows,  rlo_);
944  ClpDisjointCopyN(si->getRowUpper(), nrows,  rup_);
945  int i;
946  for (i=0;i<ncols_in;i++) 
947    originalColumn_[i]=i;
948  for (i=0;i<nrows_in;i++) 
949    originalRow_[i]=i;
950  sol_=NULL;
951  rowduals_=NULL;
952  acts_=NULL;
953
954  rcosts_=NULL;
955  colstat_=NULL;
956  rowstat_=NULL;
957}
958
959// I am not familiar enough with CoinPackedMatrix to be confident
960// that I will implement a row-ordered version of toColumnOrderedGapFree
961// properly.
962static bool isGapFree(const CoinPackedMatrix& matrix)
963{
964  const CoinBigIndex * start = matrix.getVectorStarts();
965  const int * length = matrix.getVectorLengths();
966  int i;
967  for (i = matrix.getSizeVectorLengths() - 1; i >= 0; --i) {
968    if (start[i+1] - start[i] != length[i])
969      break;
970  }
971  return (! (i >= 0));
972}
973#if     DEBUG_PRESOLVE
974static void matrix_bounds_ok(const double *lo, const double *up, int n)
975{
976  int i;
977  for (i=0; i<n; i++) {
978    PRESOLVEASSERT(lo[i] <= up[i]);
979    PRESOLVEASSERT(lo[i] < PRESOLVE_INF);
980    PRESOLVEASSERT(-PRESOLVE_INF < up[i]);
981  }
982}
983#endif
984CoinPresolveMatrix::CoinPresolveMatrix(int ncols0_in,
985                                     double maxmin_,
986                                     // end prepost members
987
988                                     ClpSimplex * si,
989
990                                     // rowrep
991                                     int nrows_in,
992                                     CoinBigIndex nelems_in,
993                               bool doStatus,
994                               double nonLinearValue) :
995
996  CoinPrePostsolveMatrix(si,
997                        ncols0_in, nrows_in, nelems_in),
998  clink_(new presolvehlink[ncols0_in+1]),
999  rlink_(new presolvehlink[nrows_in+1]),
1000
1001  dobias_(0.0),
1002
1003  nrows_(si->getNumRows()),
1004
1005  // temporary init
1006  mrstrt_(new CoinBigIndex[nrows_in+1]),
1007  hinrow_(new int[nrows_in+1]),
1008  rowels_(new double[2*nelems_in]),
1009  hcol_(new int[2*nelems_in]),
1010  integerType_(new char[ncols0_in]),
1011  feasibilityTolerance_(0.0),
1012  status_(-1),
1013  rowsToDo_(new int [nrows_in]),
1014  numberRowsToDo_(0),
1015  nextRowsToDo_(new int[nrows_in]),
1016  numberNextRowsToDo_(0),
1017  colsToDo_(new int [ncols0_in]),
1018  numberColsToDo_(0),
1019  nextColsToDo_(new int[ncols0_in]),
1020  numberNextColsToDo_(0)
1021
1022{
1023  const int bufsize = 2*nelems_in;
1024
1025  // Set up change bits etc
1026  rowChanged_ = new unsigned char[nrows_];
1027  memset(rowChanged_,0,nrows_);
1028  colChanged_ = new unsigned char[ncols_];
1029  memset(colChanged_,0,ncols_);
1030  CoinPackedMatrix * m = si->matrix();
1031
1032  // The coefficient matrix is a big hunk of stuff.
1033  // Do the copy here to try to avoid running out of memory.
1034
1035  const CoinBigIndex * start = m->getVectorStarts();
1036  const int * length = m->getVectorLengths();
1037  const int * row = m->getIndices();
1038  const double * element = m->getElements();
1039  int icol,nel=0;
1040  mcstrt_[0]=0;
1041  for (icol=0;icol<ncols_;icol++) {
1042    int j;
1043    for (j=start[icol];j<start[icol]+length[icol];j++) {
1044      hrow_[nel]=row[j];
1045      colels_[nel++]=element[j];
1046    }
1047    mcstrt_[icol+1]=nel;
1048  }
1049  assert(mcstrt_[ncols_] == nelems_);
1050  ClpDisjointCopyN(m->getVectorLengths(),ncols_,  hincol_);
1051
1052  // same thing for row rep
1053  m = new CoinPackedMatrix();
1054  m->reverseOrderedCopyOf(*si->matrix());
1055  m->removeGaps();
1056
1057
1058  ClpDisjointCopyN(m->getVectorStarts(),  nrows_,  mrstrt_);
1059  mrstrt_[nrows_] = nelems_;
1060  ClpDisjointCopyN(m->getVectorLengths(), nrows_,  hinrow_);
1061  ClpDisjointCopyN(m->getIndices(),       nelems_, hcol_);
1062  ClpDisjointCopyN(m->getElements(),      nelems_, rowels_);
1063
1064  delete m;
1065  if (si->integerInformation()) {
1066    memcpy(integerType_,si->integerInformation(),ncols_*sizeof(char));
1067  } else {
1068    ClpFillN<char>(integerType_, ncols_, 0);
1069  }
1070
1071  // Set up prohibited bits if needed
1072  if (nonLinearValue) {
1073    anyProhibited_ = true;
1074    for (icol=0;icol<ncols_;icol++) {
1075      int j;
1076      bool nonLinearColumn = false;
1077      if (cost_[icol]==nonLinearValue)
1078        nonLinearColumn=true;
1079      for (j=mcstrt_[icol];j<mcstrt_[icol+1];j++) {
1080        if (colels_[j]==nonLinearValue) {
1081          nonLinearColumn=true;
1082          setRowProhibited(hrow_[j]);
1083        }
1084      }
1085      if (nonLinearColumn)
1086        setColProhibited(icol);
1087    }
1088  } else {
1089    anyProhibited_ = false;
1090  }
1091
1092  if (doStatus) {
1093    // allow for status and solution
1094    sol_ = new double[ncols_];
1095    memcpy(sol_,si->primalColumnSolution(),ncols_*sizeof(double));;
1096    acts_ = new double [nrows_];
1097    memcpy(acts_,si->primalRowSolution(),nrows_*sizeof(double));
1098    if (!si->statusArray())
1099      si->createStatus();
1100    colstat_ = new unsigned char [nrows_+ncols_];
1101    memcpy(colstat_,si->statusArray(),
1102           (nrows_+ncols_)*sizeof(unsigned char));
1103    rowstat_ = colstat_+ncols_;
1104  }
1105
1106  // the original model's fields are now unneeded - free them
1107 
1108  si->resize(0,0);
1109
1110#if     DEBUG_PRESOLVE
1111  matrix_bounds_ok(rlo_, rup_, nrows_);
1112  matrix_bounds_ok(clo_, cup_, ncols_);
1113#endif
1114
1115#if 0
1116  for (i=0; i<nrows; ++i)
1117    printf("NR: %6d\n", hinrow[i]);
1118  for (int i=0; i<ncols; ++i)
1119    printf("NC: %6d\n", hincol[i]);
1120#endif
1121
1122  presolve_make_memlists(mcstrt_, hincol_, clink_, ncols_);
1123  presolve_make_memlists(mrstrt_, hinrow_, rlink_, nrows_);
1124
1125  // this allows last col/row to expand up to bufsize-1 (22);
1126  // this must come after the calls to presolve_prefix
1127  mcstrt_[ncols_] = bufsize-1;
1128  mrstrt_[nrows_] = bufsize-1;
1129
1130#if     CHECK_CONSISTENCY
1131  consistent(false);
1132#endif
1133}
1134
1135void CoinPresolveMatrix::update_model(ClpSimplex * si,
1136                                     int nrows0,
1137                                     int ncols0,
1138                                     CoinBigIndex nelems0)
1139{
1140  si->loadProblem(ncols_, nrows_, mcstrt_, hrow_, colels_, hincol_,
1141                 clo_, cup_, cost_, rlo_, rup_);
1142
1143  delete [] si->integerInformation();
1144  int numberIntegers=0;
1145  for (int i=0; i<ncols_; i++) {
1146    if (integerType_[i])
1147      numberIntegers++;
1148  }
1149  if (numberIntegers) 
1150    si->copyInIntegerInformation(integerType_);
1151  else
1152    si->copyInIntegerInformation(NULL);
1153
1154#if     PRESOLVE_SUMMARY
1155  printf("NEW NCOL/NROW/NELS:  %d(-%d) %d(-%d) %d(-%d)\n",
1156         ncols_, ncols0-ncols_,
1157         nrows_, nrows0-nrows_,
1158         si->getNumElements(), nelems0-si->getNumElements());
1159#endif
1160  si->setDblParam(ClpObjOffset,originalOffset_-dobias_);
1161
1162}
1163
1164
1165
1166
1167
1168
1169
1170
1171
1172
1173
1174////////////////  POSTSOLVE
1175
1176CoinPostsolveMatrix::CoinPostsolveMatrix(ClpSimplex*  si,
1177                                       int ncols0_in,
1178                                       int nrows0_in,
1179                                       CoinBigIndex nelems0,
1180                                   
1181                                       double maxmin_,
1182                                       // end prepost members
1183
1184                                       double *sol_in,
1185                                       double *acts_in,
1186
1187                                       unsigned char *colstat_in,
1188                                       unsigned char *rowstat_in) :
1189  CoinPrePostsolveMatrix(si,
1190                        ncols0_in, nrows0_in, nelems0),
1191
1192  free_list_(0),
1193  // link, free_list, maxlink
1194  maxlink_(2*nelems0),
1195  link_(new int[/*maxlink*/ 2*nelems0]),
1196     
1197  cdone_(new char[ncols0_]),
1198  rdone_(new char[nrows0_in]),
1199
1200  nrows_(si->getNumRows()),
1201  nrows0_(nrows0_in)
1202{
1203
1204  sol_=sol_in;
1205  rowduals_=NULL;
1206  acts_=acts_in;
1207
1208  rcosts_=NULL;
1209  colstat_=colstat_in;
1210  rowstat_=rowstat_in;
1211
1212  // this is the *reduced* model, which is probably smaller
1213  int ncols1 = si->getNumCols();
1214  int nrows1 = si->getNumRows();
1215
1216  const CoinPackedMatrix * m = si->matrix();
1217
1218  if (! isGapFree(*m)) {
1219    CoinPresolveAction::throwCoinError("Matrix not gap free",
1220                                      "CoinPostsolveMatrix");
1221  }
1222
1223  const CoinBigIndex nelemsr = m->getNumElements();
1224
1225  ClpDisjointCopyN(m->getVectorStarts(), ncols1, mcstrt_);
1226  mcstrt_[ncols_] = nelems0;    // ??
1227  ClpDisjointCopyN(m->getVectorLengths(),ncols1,  hincol_);
1228  ClpDisjointCopyN(m->getIndices(),      nelemsr, hrow_);
1229  ClpDisjointCopyN(m->getElements(),     nelemsr, colels_);
1230
1231
1232#if     0 && DEBUG_PRESOLVE
1233  presolve_check_costs(model, &colcopy);
1234#endif
1235
1236  // This determines the size of the data structure that contains
1237  // the matrix being postsolved.  Links are taken from the free_list
1238  // to recreate matrix entries that were presolved away,
1239  // and links are added to the free_list when entries created during
1240  // presolve are discarded.  There is never a need to gc this list.
1241  // Naturally, it should contain
1242  // exactly nelems0 entries "in use" when postsolving is done,
1243  // but I don't know whether the matrix could temporarily get
1244  // larger during postsolving.  Substitution into more than two
1245  // rows could do that, in principle.  I am being very conservative
1246  // here by reserving much more than the amount of space I probably need.
1247  // If this guess is wrong, check_free_list may be called.
1248  //  int bufsize = 2*nelems0;
1249
1250  memset(cdone_, -1, ncols0_);
1251  memset(rdone_, -1, nrows0_);
1252
1253  rowduals_ = new double[nrows0_];
1254  ClpDisjointCopyN(si->getRowPrice(), nrows1, rowduals_);
1255
1256  rcosts_ = new double[ncols0_];
1257  ClpDisjointCopyN(si->getReducedCost(), ncols1, rcosts_);
1258  if (maxmin_<0.0) {
1259    // change so will look as if minimize
1260    int i;
1261    for (i=0;i<nrows1;i++)
1262      rowduals_[i] = - rowduals_[i];
1263    for (i=0;i<ncols1;i++) {
1264      rcosts_[i] = - rcosts_[i];
1265    }
1266  }
1267
1268  //ClpDisjointCopyN(si->getRowUpper(), nrows1, rup_);
1269  //ClpDisjointCopyN(si->getRowLower(), nrows1, rlo_);
1270
1271  ClpDisjointCopyN(si->getColSolution(), ncols1, sol_);
1272  si->setDblParam(ClpObjOffset,originalOffset_);
1273
1274  for (int j=0; j<ncols1; j++) {
1275    CoinBigIndex kcs = mcstrt_[j];
1276    CoinBigIndex kce = kcs + hincol_[j];
1277    for (CoinBigIndex k=kcs; k<kce; ++k) {
1278      link_[k] = k+1;
1279    }
1280  }
1281  {
1282    int ml = maxlink_;
1283    for (CoinBigIndex k=nelemsr; k<ml; ++k)
1284      link_[k] = k+1;
1285    link_[ml-1] = NO_LINK;
1286  }
1287  free_list_ = nelemsr;
1288}
1289/* This is main part of Presolve */
1290ClpSimplex * 
1291ClpPresolve::gutsOfPresolvedModel(ClpSimplex * originalModel,
1292                                  double feasibilityTolerance,
1293                                  bool keepIntegers,
1294                                  int numberPasses,
1295                                  bool dropNames)
1296{
1297  ncols_ = originalModel->getNumCols();
1298  nrows_ = originalModel->getNumRows();
1299  nelems_ = originalModel->getNumElements();
1300  numberPasses_ = numberPasses;
1301
1302  double maxmin = originalModel->getObjSense();
1303  originalModel_ = originalModel;
1304  delete [] originalColumn_;
1305  originalColumn_ = new int[ncols_];
1306  delete [] originalRow_;
1307  originalRow_ = new int[nrows_];
1308
1309  // result is 0 - okay, 1 infeasible, -1 go round again
1310  int result = -1;
1311 
1312  // User may have deleted - its their responsibility
1313  presolvedModel_=NULL;
1314  // Messages
1315  CoinMessages messages = CoinMessage(originalModel->messages().language());
1316  while (result==-1) {
1317
1318    // make new copy
1319    if (saveFile_=="") {
1320      delete presolvedModel_;
1321      presolvedModel_ = new ClpSimplex(*originalModel);
1322    } else {
1323      presolvedModel_=originalModel;
1324    }
1325    presolvedModel_->dropNames();
1326
1327    // drop integer information if wanted
1328    if (!keepIntegers)
1329      presolvedModel_->deleteIntegerInformation();
1330
1331   
1332    CoinPresolveMatrix prob(ncols_,
1333                        maxmin,
1334                        presolvedModel_,
1335                        nrows_, nelems_,true,nonLinearValue_);
1336    // make sure row solution correct
1337    {
1338      double *colels    = prob.colels_;
1339      int *hrow         = prob.hrow_;
1340      CoinBigIndex *mcstrt              = prob.mcstrt_;
1341      int *hincol               = prob.hincol_;
1342      int ncols         = prob.ncols_;
1343     
1344     
1345      double * csol = prob.sol_;
1346      double * acts = prob.acts_;
1347      int nrows = prob.nrows_;
1348
1349      int colx;
1350
1351      memset(acts,0,nrows*sizeof(double));
1352     
1353      for (colx = 0; colx < ncols; ++colx) {
1354        double solutionValue = csol[colx];
1355        for (int i=mcstrt[colx]; i<mcstrt[colx]+hincol[colx]; ++i) {
1356          int row = hrow[i];
1357          double coeff = colels[i];
1358          acts[row] += solutionValue*coeff;
1359        }
1360      }
1361    }
1362    prob.handler_ = presolvedModel_->messageHandler();
1363    prob.messages_ = CoinMessage(presolvedModel_->messages().language());
1364
1365    // move across feasibility tolerance
1366    prob.feasibilityTolerance_ = feasibilityTolerance;
1367
1368    // Do presolve
1369
1370    paction_ = presolve(&prob);
1371
1372    result =0; 
1373
1374    if (prob.status_==0&&paction_) {
1375      // Looks feasible but double check to see if anything slipped through
1376      int n             = prob.ncols_;
1377      double * lo = prob.clo_;
1378      double * up = prob.cup_;
1379      int i;
1380     
1381      for (i=0;i<n;i++) {
1382        if (up[i]<lo[i]) {
1383          if (up[i]<lo[i]-1.0e-8) {
1384            // infeasible
1385            prob.status_=1;
1386          } else {
1387            up[i]=lo[i];
1388          }
1389        }
1390      }
1391     
1392      n = prob.nrows_;
1393      lo = prob.rlo_;
1394      up = prob.rup_;
1395
1396      for (i=0;i<n;i++) {
1397        if (up[i]<lo[i]) {
1398          if (up[i]<lo[i]-1.0e-8) {
1399            // infeasible
1400            prob.status_=1;
1401          } else {
1402            up[i]=lo[i];
1403          }
1404        }
1405      }
1406    }
1407    if (prob.status_==0&&paction_) {
1408      // feasible
1409   
1410      prob.update_model(presolvedModel_, nrows_, ncols_, nelems_);
1411      // copy status and solution
1412      memcpy(presolvedModel_->primalColumnSolution(),
1413             prob.sol_,prob.ncols_*sizeof(double));
1414      memcpy(presolvedModel_->primalRowSolution(),
1415             prob.acts_,prob.nrows_*sizeof(double));
1416      memcpy(presolvedModel_->statusArray(),
1417             prob.colstat_,prob.ncols_*sizeof(unsigned char));
1418      memcpy(presolvedModel_->statusArray()+prob.ncols_,
1419             prob.rowstat_,prob.nrows_*sizeof(unsigned char));
1420      delete [] prob.sol_;
1421      delete [] prob.acts_;
1422      delete [] prob.colstat_;
1423      prob.sol_=NULL;
1424      prob.acts_=NULL;
1425      prob.colstat_=NULL;
1426     
1427      int ncolsNow = presolvedModel_->getNumCols();
1428      memcpy(originalColumn_,prob.originalColumn_,ncolsNow*sizeof(int));
1429      delete [] prob.originalColumn_;
1430      prob.originalColumn_=NULL;
1431      int nrowsNow = presolvedModel_->getNumRows();
1432      memcpy(originalRow_,prob.originalRow_,nrowsNow*sizeof(int));
1433      delete [] prob.originalRow_;
1434      prob.originalRow_=NULL;
1435      if (!dropNames&&originalModel->lengthNames()) {
1436        // Redo names
1437        int iRow;
1438        std::vector<std::string> rowNames;
1439        rowNames.reserve(nrowsNow);
1440        for (iRow=0;iRow<nrowsNow;iRow++) {
1441          int kRow = originalRow_[iRow];
1442          rowNames.push_back(originalModel->rowName(kRow));
1443        }
1444     
1445        int iColumn;
1446        std::vector<std::string> columnNames;
1447        columnNames.reserve(ncolsNow);
1448        for (iColumn=0;iColumn<ncolsNow;iColumn++) {
1449          int kColumn = originalColumn_[iColumn];
1450          columnNames.push_back(originalModel->columnName(kColumn));
1451        }
1452        presolvedModel_->copyNames(rowNames,columnNames);
1453      }
1454      // now clean up integer variables.  This can modify original
1455      int i;
1456      const char * information = presolvedModel_->integerInformation();
1457      if (information) {
1458        int numberChanges=0;
1459        double * lower0 = originalModel_->columnLower();
1460        double * upper0 = originalModel_->columnUpper();
1461        double * lower = presolvedModel_->columnLower();
1462        double * upper = presolvedModel_->columnUpper();
1463        for (i=0;i<ncolsNow;i++) {
1464          if (!information[i])
1465            continue;
1466          int iOriginal = originalColumn_[i];
1467          double lowerValue0 = lower0[iOriginal];
1468          double upperValue0 = upper0[iOriginal];
1469          double lowerValue = ceil(lower[i]-1.0e-5);
1470          double upperValue = floor(upper[i]+1.0e-5);
1471          lower[i]=lowerValue;
1472          upper[i]=upperValue;
1473          if (lowerValue>upperValue) {
1474            numberChanges++;
1475            presolvedModel_->messageHandler()->message(COIN_PRESOLVE_COLINFEAS,
1476                                                       messages)
1477                                                         <<iOriginal
1478                                                         <<lowerValue
1479                                                         <<upperValue
1480                                                         <<CoinMessageEol;
1481            result=1;
1482          } else {
1483            if (lowerValue>lowerValue0+1.0e-8) {
1484              lower0[iOriginal] = lowerValue;
1485              numberChanges++;
1486            }
1487            if (upperValue<upperValue0-1.0e-8) {
1488              upper0[iOriginal] = upperValue;
1489              numberChanges++;
1490            }
1491          }       
1492        }
1493        if (numberChanges) {
1494          presolvedModel_->messageHandler()->message(COIN_PRESOLVE_INTEGERMODS,
1495                                                     messages)
1496                                                       <<numberChanges
1497                                                       <<CoinMessageEol;
1498          if (!result) {
1499            result = -1; // round again
1500          }
1501        }
1502      }
1503    } else {
1504      // infeasible
1505      delete [] prob.sol_;
1506      delete [] prob.acts_;
1507      delete [] prob.colstat_;
1508      result=1;
1509    }
1510  }
1511  if (!result) {
1512    int nrowsAfter = presolvedModel_->getNumRows();
1513    int ncolsAfter = presolvedModel_->getNumCols();
1514    CoinBigIndex nelsAfter = presolvedModel_->getNumElements();
1515    presolvedModel_->messageHandler()->message(COIN_PRESOLVE_STATS,
1516                                               messages)
1517                                                 <<nrowsAfter<< -(nrows_ - nrowsAfter)
1518                                                 << ncolsAfter<< -(ncols_ - ncolsAfter)
1519                                                 <<nelsAfter<< -(nelems_ - nelsAfter)
1520                                                 <<CoinMessageEol;
1521  } else {
1522    gutsOfDestroy();
1523    delete presolvedModel_;
1524    presolvedModel_=NULL;
1525  }
1526  return presolvedModel_;
1527}
1528
1529
Note: See TracBrowser for help on using the repository browser.