source: trunk/Clp/src/ClpPresolve.cpp @ 1370

Last change on this file since 1370 was 1370, checked in by forrest, 10 years ago

add ids

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