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

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

trying to make faster

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