source: trunk/ClpPresolve.cpp @ 650

Last change on this file since 650 was 650, checked in by forrest, 14 years ago

stuff

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