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

Last change on this file since 1344 was 1344, checked in by forrest, 11 years ago

changes to simplex and lots of stuff and start Mumps cholesky

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