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

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

out compiler warnings and stability improvements

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 49.2 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#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      //#define IMPLIED 3
509#ifdef IMPLIED
510      int fill_level=3;
511#define IMPLIED2 99
512#if IMPLIED!=3
513#if IMPLIED>2&&IMPLIED<11
514      fill_level=IMPLIED;
515      printf("** fill_level == %d !\n",fill_level);
516#endif
517#if IMPLIED>11&&IMPLIED<21
518      fill_level=-(IMPLIED-10);
519      printf("** fill_level == %d !\n",fill_level);
520#endif
521#endif
522#else
523      int fill_level=2;
524#endif
525      int whichPass=0;
526      while (1) {
527        whichPass++;
528        prob->pass_++;
529        const CoinPresolveAction * const paction1 = paction_;
530
531        if (slackd) {
532          bool notFinished = true;
533          while (notFinished) 
534            paction_ = slack_doubleton_action::presolve(prob, paction_,
535                                                        notFinished);
536          if (prob->status_)
537            break;
538        }
539        if (dual&&whichPass==1) {
540          // this can also make E rows so do one bit here
541          paction_ = remove_dual_action::presolve(prob, paction_);
542          if (prob->status_)
543            break;
544        }
545
546        if (doubleton) {
547          paction_ = doubleton_action::presolve(prob, paction_);
548          if (prob->status_)
549            break;
550        }
551        if (tripleton) {
552          paction_ = tripleton_action::presolve(prob, paction_);
553          if (prob->status_)
554            break;
555        }
556
557        if (zerocost) {
558          paction_ = do_tighten_action::presolve(prob, paction_);
559          if (prob->status_)
560            break;
561        }
562#ifndef NO_FORCING
563        if (forcing) {
564          paction_ = forcing_constraint_action::presolve(prob, paction_);
565          if (prob->status_)
566            break;
567        }
568#endif
569
570        if (ifree&&(whichPass%5)==1) {
571          paction_ = implied_free_action::presolve(prob, paction_,fill_level);
572        if (prob->status_)
573          break;
574        }
575
576#if     PRESOLVE_DEBUG
577        check_sol(prob,1.0e0);
578#endif
579
580#if     PRESOLVE_CONSISTENCY
581//      presolve_links_ok(prob->rlink_, prob->mrstrt_, prob->hinrow_,
582//                        prob->nrows_);
583        presolve_links_ok(prob,false,true) ;
584#endif
585
586//#if   PRESOLVE_DEBUG
587//      presolve_no_zeros(prob->mcstrt_, prob->colels_, prob->hincol_,
588//                        prob->ncols_);
589//#endif
590//#if   PRESOLVE_CONSISTENCY
591//      prob->consistent();
592//#endif
593#if     PRESOLVE_CONSISTENCY
594        presolve_no_zeros(prob,true,false) ;
595        presolve_consistent(prob,true) ;
596#endif
597
598         
599        // set up for next pass
600        // later do faster if many changes i.e. memset and memcpy
601        prob->numberRowsToDo_ = prob->numberNextRowsToDo_;
602        int kcheck;
603        bool found=false;
604        kcheck=-1;
605        for (i=0;i<prob->numberNextRowsToDo_;i++) {
606          int index = prob->nextRowsToDo_[i];
607          prob->unsetRowChanged(index);
608          prob->rowsToDo_[i] = index;
609          if (index==kcheck) {
610            printf("row %d on list after pass %d\n",kcheck,
611                   whichPass);
612            found=true;
613          }
614        }
615        if (!found&&kcheck>=0)
616          prob->rowsToDo_[prob->numberRowsToDo_++]=kcheck;
617        prob->numberNextRowsToDo_=0;
618        prob->numberColsToDo_ = prob->numberNextColsToDo_;
619        kcheck=-1;
620        found=false;
621        for (i=0;i<prob->numberNextColsToDo_;i++) {
622          int index = prob->nextColsToDo_[i];
623          prob->unsetColChanged(index);
624          prob->colsToDo_[i] = index;
625          if (index==kcheck) {
626            printf("col %d on list after pass %d\n",kcheck,
627                   whichPass);
628            found=true;
629          }
630        }
631        if (!found&&kcheck>=0)
632          prob->colsToDo_[prob->numberColsToDo_++]=kcheck;
633        prob->numberNextColsToDo_=0;
634        if (paction_ == paction1&&fill_level>0)
635          break;
636      }
637      // say look at all
638      int i;
639      if (!prob->anyProhibited()) {
640        for (i=0;i<nrows_;i++) 
641          prob->rowsToDo_[i]=i;
642        prob->numberRowsToDo_=nrows_;
643        for (i=0;i<ncols_;i++) 
644          prob->colsToDo_[i]=i;
645        prob->numberColsToDo_=ncols_;
646      } else {
647        // some stuff must be left alone
648        prob->numberRowsToDo_=0;
649        for (i=0;i<nrows_;i++) 
650          if (!prob->rowProhibited(i))
651            prob->rowsToDo_[prob->numberRowsToDo_++]=i;
652        prob->numberColsToDo_=0;
653        for (i=0;i<ncols_;i++) 
654          if (!prob->colProhibited(i))
655            prob->colsToDo_[prob->numberColsToDo_++]=i;
656      }
657      // now expensive things
658      // this caused world.mps to run into numerical difficulties
659#ifdef PRESOLVE_SUMMARY
660      printf("Starting expensive\n");
661#endif
662
663      if (dual) {
664        int itry;
665        for (itry=0;itry<5;itry++) {
666          paction_ = remove_dual_action::presolve(prob, paction_);
667          if (prob->status_)
668            break;
669          const CoinPresolveAction * const paction2 = paction_;
670          if (ifree) {
671#ifdef IMPLIED
672#if IMPLIED2 ==0
673            int fill_level=0; // switches off substitution
674#elif IMPLIED2!=99
675            int fill_level=IMPLIED2;
676#endif
677#endif
678            if ((itry&1)==0)
679              paction_ = implied_free_action::presolve(prob, paction_,fill_level);
680            if (prob->status_)
681              break;
682          }
683          if (paction_ == paction2)
684            break;
685        }
686      } else if (ifree) {
687        // just free
688#ifdef IMPLIED
689#if IMPLIED2 ==0
690        int fill_level=0; // switches off substitution
691#elif IMPLIED2!=99
692        int fill_level=IMPLIED2;
693#endif
694#endif
695        paction_ = implied_free_action::presolve(prob, paction_,fill_level);
696        if (prob->status_)
697          break;
698      }
699#if     PRESOLVE_DEBUG
700      check_sol(prob,1.0e0);
701#endif
702      if (dupcol) {
703        // maybe allow integer columns to be checked
704        if ((presolveActions_&512)!=0)
705          prob->setPresolveOptions(prob->presolveOptions()|1);
706        paction_ = dupcol_action::presolve(prob, paction_);
707        if (prob->status_)
708          break;
709      }
710#if     PRESOLVE_DEBUG
711        check_sol(prob,1.0e0);
712#endif
713     
714      if (duprow) {
715        paction_ = duprow_action::presolve(prob, paction_);
716        if (prob->status_)
717          break;
718      }
719#if     PRESOLVE_DEBUG
720      check_sol(prob,1.0e0);
721#endif
722      bool stopLoop=false;
723      {
724        int * hinrow = prob->hinrow_;
725        int numberDropped=0;
726        for (i=0;i<nrows_;i++) 
727          if (!hinrow[i])
728            numberDropped++;
729
730        prob->messageHandler()->message(COIN_PRESOLVE_PASS,
731                                        messages)
732                                          <<numberDropped<<iLoop+1
733                                          <<CoinMessageEol;
734        //printf("%d rows dropped after pass %d\n",numberDropped,
735        //     iLoop+1);
736        if (numberDropped==lastDropped)
737          stopLoop=true;
738        else
739          lastDropped = numberDropped;
740      }
741      // Do this here as not very loopy
742      if (slackSingleton) {
743        // On most passes do not touch costed slacks
744        if (paction_ != paction0&&!stopLoop) {
745          paction_ = slack_singleton_action::presolve(prob, paction_,NULL);
746        } else {
747          // do costed if Clp (at end as ruins rest of presolve)
748          paction_ = slack_singleton_action::presolve(prob, paction_,rowObjective_);
749          stopLoop=true;
750        }
751      }
752#if     PRESOLVE_DEBUG
753      check_sol(prob,1.0e0);
754#endif
755      if (paction_ == paction0||stopLoop)
756        break;
757         
758    }
759  }
760  if (!prob->status_) {
761    paction_ = drop_zero_coefficients(prob, paction_);
762#if     PRESOLVE_DEBUG
763        check_sol(prob,1.0e0);
764#endif
765
766    paction_ = drop_empty_cols_action::presolve(prob, paction_);
767    paction_ = drop_empty_rows_action::presolve(prob, paction_);
768#if     PRESOLVE_DEBUG
769        check_sol(prob,1.0e0);
770#endif
771  }
772 
773  if (prob->status_) {
774    if (prob->status_==1) 
775          prob->messageHandler()->message(COIN_PRESOLVE_INFEAS,
776                                             messages)
777                                               <<prob->feasibilityTolerance_
778                                               <<CoinMessageEol;
779    else if (prob->status_==2)
780          prob->messageHandler()->message(COIN_PRESOLVE_UNBOUND,
781                                             messages)
782                                               <<CoinMessageEol;
783    else
784          prob->messageHandler()->message(COIN_PRESOLVE_INFEASUNBOUND,
785                                             messages)
786                                               <<CoinMessageEol;
787    // get rid of data
788    destroyPresolve();
789  }
790  return (paction_);
791}
792
793void check_djs(CoinPostsolveMatrix *prob);
794
795
796// We could have implemented this by having each postsolve routine
797// directly call the next one, but this may make it easier to add debugging checks.
798void ClpPresolve::postsolve(CoinPostsolveMatrix &prob)
799{
800  {
801    // Check activities
802    double *colels      = prob.colels_;
803    int *hrow           = prob.hrow_;
804    CoinBigIndex *mcstrt                = prob.mcstrt_;
805    int *hincol         = prob.hincol_;
806    int *link           = prob.link_;
807    int ncols           = prob.ncols_;
808
809    char *cdone = prob.cdone_;
810
811    double * csol = prob.sol_;
812    int nrows = prob.nrows_;
813
814    int colx;
815
816    double * rsol = prob.acts_;
817    memset(rsol,0,nrows*sizeof(double));
818
819    for (colx = 0; colx < ncols; ++colx) {
820      if (cdone[colx]) {
821        CoinBigIndex k = mcstrt[colx];
822        int nx = hincol[colx];
823        double solutionValue = csol[colx];
824        for (int i=0; i<nx; ++i) {
825          int row = hrow[k];
826          double coeff = colels[k];
827          k = link[k];
828          rsol[row] += solutionValue*coeff;
829        }
830      }
831    }
832  }
833  const CoinPresolveAction *paction = paction_;
834  //#define PRESOLVE_DEBUG 1
835#if     PRESOLVE_DEBUG
836  // Check only works after first one
837  int checkit=-1;
838#endif
839 
840  while (paction) {
841#if PRESOLVE_DEBUG
842    printf("POSTSOLVING %s\n", paction->name());
843#endif
844
845    paction->postsolve(&prob);
846   
847#if     PRESOLVE_DEBUG
848  {
849    int nr=0;
850    int i;
851    for (i=0;i<prob.nrows_;i++) {
852      if ((prob.rowstat_[i]&7)==1)
853        nr++;
854    }
855    int nc=0;
856    for (i=0;i<prob.ncols_;i++) {
857      if ((prob.colstat_[i]&7)==1)
858        nc++;
859    }
860    printf("%d rows (%d basic), %d cols (%d basic)\n",prob.nrows_,nr,prob.ncols_,nc);
861  }
862    checkit++;
863    if (prob.colstat_&&checkit>0) {
864      presolve_check_nbasic(&prob) ;
865      presolve_check_sol(&prob,2,2,1) ;
866    }
867#endif
868    paction = paction->next;
869#if     PRESOLVE_DEBUG
870//  check_djs(&prob);
871    if (checkit>0)
872      presolve_check_reduced_costs(&prob) ;
873#endif
874  }   
875#if     PRESOLVE_DEBUG
876  if (prob.colstat_) {
877    presolve_check_nbasic(&prob) ;
878    presolve_check_sol(&prob,2,2,1) ;
879  }
880#endif
881#undef PRESOLVE_DEBUG
882 
883#if     0 && PRESOLVE_DEBUG
884  for (i=0; i<ncols0; i++) {
885    if (!cdone[i]) {
886      printf("!cdone[%d]\n", i);
887      abort();
888    }
889  }
890 
891  for (i=0; i<nrows0; i++) {
892    if (!rdone[i]) {
893      printf("!rdone[%d]\n", i);
894      abort();
895    }
896  }
897 
898 
899  for (i=0; i<ncols0; i++) {
900    if (sol[i] < -1e10 || sol[i] > 1e10)
901      printf("!!!%d %g\n", i, sol[i]);
902   
903  }
904 
905 
906#endif
907 
908#if     0 && PRESOLVE_DEBUG
909  // debug check:  make sure we ended up with same original matrix
910  {
911    int identical = 1;
912   
913    for (int i=0; i<ncols0; i++) {
914      PRESOLVEASSERT(hincol[i] == &prob->mcstrt0[i+1] - &prob->mcstrt0[i]);
915      CoinBigIndex kcs0 = &prob->mcstrt0[i];
916      CoinBigIndex kcs = mcstrt[i];
917      int n = hincol[i];
918      for (int k=0; k<n; k++) {
919        CoinBigIndex k1 = presolve_find_row1(&prob->hrow0[kcs0+k], kcs, kcs+n, hrow);
920
921        if (k1 == kcs+n) {
922          printf("ROW %d NOT IN COL %d\n", &prob->hrow0[kcs0+k], i);
923          abort();
924        }
925
926        if (colels[k1] != &prob->dels0[kcs0+k])
927          printf("BAD COLEL[%d %d %d]:  %g\n",
928                 k1, i, &prob->hrow0[kcs0+k], colels[k1] - &prob->dels0[kcs0+k]);
929
930        if (kcs0+k != k1)
931          identical=0;
932      }
933    }
934    printf("identical? %d\n", identical);
935  }
936#endif
937}
938
939
940
941
942
943
944
945
946static inline double getTolerance(const ClpSimplex  *si, ClpDblParam key)
947{
948  double tol;
949  if (! si->getDblParam(key, tol)) {
950    CoinPresolveAction::throwCoinError("getDblParam failed",
951                                      "CoinPrePostsolveMatrix::CoinPrePostsolveMatrix");
952  }
953  return (tol);
954}
955
956
957// Assumptions:
958// 1. nrows>=m.getNumRows()
959// 2. ncols>=m.getNumCols()
960//
961// In presolve, these values are equal.
962// In postsolve, they may be inequal, since the reduced problem
963// may be smaller, but we need room for the large problem.
964// ncols may be larger than si.getNumCols() in postsolve,
965// this at that point si will be the reduced problem,
966// but we need to reserve enough space for the original problem.
967CoinPrePostsolveMatrix::CoinPrePostsolveMatrix(const ClpSimplex * si,
968                                             int ncols_in,
969                                             int nrows_in,
970                                             CoinBigIndex nelems_in,
971                                               double bulkRatio) 
972  : ncols_(si->getNumCols()),
973    nrows_(si->getNumRows()),
974    nelems_(si->getNumElements()),
975    ncols0_(ncols_in),
976    nrows0_(nrows_in),
977    bulkRatio_(bulkRatio),
978    mcstrt_(new CoinBigIndex[ncols_in+1]),
979    hincol_(new int[ncols_in+1]),
980    cost_(new double[ncols_in]),
981    clo_(new double[ncols_in]),
982    cup_(new double[ncols_in]),
983    rlo_(new double[nrows_in]),
984    rup_(new double[nrows_in]),
985    originalColumn_(new int[ncols_in]),
986    originalRow_(new int[nrows_in]),
987    ztolzb_(getTolerance(si, ClpPrimalTolerance)),
988    ztoldj_(getTolerance(si, ClpDualTolerance)),
989    maxmin_(si->getObjSense()),
990    sol_(NULL),
991    rowduals_(NULL),
992    acts_(NULL),
993    rcosts_(NULL),
994    colstat_(NULL),
995    rowstat_(NULL),
996    handler_(NULL),
997    defaultHandler_(false),
998    messages_()
999
1000{
1001  bulk0_ = static_cast<CoinBigIndex> (bulkRatio_*nelems_in);
1002  hrow_  = new int   [bulk0_];
1003  colels_ = new double[bulk0_];
1004  si->getDblParam(ClpObjOffset,originalOffset_);
1005  int ncols = si->getNumCols();
1006  int nrows = si->getNumRows();
1007
1008  setMessageHandler(si->messageHandler()) ;
1009
1010  ClpDisjointCopyN(si->getColLower(), ncols, clo_);
1011  ClpDisjointCopyN(si->getColUpper(), ncols, cup_);
1012  //ClpDisjointCopyN(si->getObjCoefficients(), ncols, cost_);
1013  double offset;
1014  ClpDisjointCopyN(si->objectiveAsObject()->gradient(si,si->getColSolution(),offset,true), ncols, cost_);
1015  ClpDisjointCopyN(si->getRowLower(), nrows,  rlo_);
1016  ClpDisjointCopyN(si->getRowUpper(), nrows,  rup_);
1017  int i;
1018  for (i=0;i<ncols_in;i++) 
1019    originalColumn_[i]=i;
1020  for (i=0;i<nrows_in;i++) 
1021    originalRow_[i]=i;
1022  sol_=NULL;
1023  rowduals_=NULL;
1024  acts_=NULL;
1025
1026  rcosts_=NULL;
1027  colstat_=NULL;
1028  rowstat_=NULL;
1029}
1030
1031// I am not familiar enough with CoinPackedMatrix to be confident
1032// that I will implement a row-ordered version of toColumnOrderedGapFree
1033// properly.
1034static bool isGapFree(const CoinPackedMatrix& matrix)
1035{
1036  const CoinBigIndex * start = matrix.getVectorStarts();
1037  const int * length = matrix.getVectorLengths();
1038  int i = matrix.getSizeVectorLengths() - 1;
1039  // Quick check
1040  if (matrix.getNumElements()==start[i]) {
1041    return true;
1042  } else {
1043    for (i = matrix.getSizeVectorLengths() - 1; i >= 0; --i) {
1044      if (start[i+1] - start[i] != length[i])
1045        break;
1046    }
1047    return (! (i >= 0));
1048  }
1049}
1050#if     PRESOLVE_DEBUG
1051static void matrix_bounds_ok(const double *lo, const double *up, int n)
1052{
1053  int i;
1054  for (i=0; i<n; i++) {
1055    PRESOLVEASSERT(lo[i] <= up[i]);
1056    PRESOLVEASSERT(lo[i] < PRESOLVE_INF);
1057    PRESOLVEASSERT(-PRESOLVE_INF < up[i]);
1058  }
1059}
1060#endif
1061CoinPresolveMatrix::CoinPresolveMatrix(int ncols0_in,
1062                                     double maxmin,
1063                                     // end prepost members
1064
1065                                     ClpSimplex * si,
1066
1067                                     // rowrep
1068                                     int nrows_in,
1069                                     CoinBigIndex nelems_in,
1070                                       bool doStatus,
1071                                       double nonLinearValue,
1072                                       double bulkRatio) :
1073
1074  CoinPrePostsolveMatrix(si,
1075                        ncols0_in, nrows_in, nelems_in,bulkRatio),
1076  clink_(new presolvehlink[ncols0_in+1]),
1077  rlink_(new presolvehlink[nrows_in+1]),
1078
1079  dobias_(0.0),
1080
1081
1082  // temporary init
1083  integerType_(new unsigned char[ncols0_in]),
1084  tuning_(false),
1085  startTime_(0.0),
1086  feasibilityTolerance_(0.0),
1087  status_(-1),
1088  colsToDo_(new int [ncols0_in]),
1089  numberColsToDo_(0),
1090  nextColsToDo_(new int[ncols0_in]),
1091  numberNextColsToDo_(0),
1092  rowsToDo_(new int [nrows_in]),
1093  numberRowsToDo_(0),
1094  nextRowsToDo_(new int[nrows_in]),
1095  numberNextRowsToDo_(0),
1096  presolveOptions_(0)
1097{
1098  const int bufsize = bulk0_;
1099
1100  nrows_ = si->getNumRows() ;
1101
1102  // Set up change bits etc
1103  rowChanged_ = new unsigned char[nrows_];
1104  memset(rowChanged_,0,nrows_);
1105  colChanged_ = new unsigned char[ncols_];
1106  memset(colChanged_,0,ncols_);
1107  CoinPackedMatrix * m = si->matrix();
1108
1109  // The coefficient matrix is a big hunk of stuff.
1110  // Do the copy here to try to avoid running out of memory.
1111
1112  const CoinBigIndex * start = m->getVectorStarts();
1113  const int * row = m->getIndices();
1114  const double * element = m->getElements();
1115  int icol,nel=0;
1116  mcstrt_[0]=0;
1117  ClpDisjointCopyN(m->getVectorLengths(),ncols_,  hincol_);
1118  for (icol=0;icol<ncols_;icol++) {
1119    CoinBigIndex j;
1120    for (j=start[icol];j<start[icol]+hincol_[icol];j++) {
1121      hrow_[nel]=row[j];
1122      if (fabs(element[j])>ZTOLDP)
1123        colels_[nel++]=element[j];
1124    }
1125    mcstrt_[icol+1]=nel;
1126    hincol_[icol]=nel-mcstrt_[icol];
1127  }
1128
1129  // same thing for row rep
1130  CoinPackedMatrix * mRow = new CoinPackedMatrix();
1131  mRow->setExtraGap(0.0);
1132  mRow->setExtraMajor(0.0);
1133  mRow->reverseOrderedCopyOf(*m);
1134  //mRow->removeGaps();
1135  //mRow->setExtraGap(0.0);
1136
1137  // Now get rid of matrix
1138  si->createEmptyMatrix();
1139
1140  double * el = mRow->getMutableElements();
1141  int * ind = mRow->getMutableIndices();
1142  CoinBigIndex * strt = mRow->getMutableVectorStarts();
1143  int * len = mRow->getMutableVectorLengths();
1144  // Do carefully to save memory
1145  rowels_ = new double[bulk0_];
1146  ClpDisjointCopyN(el,      nelems_, rowels_);
1147  mRow->nullElementArray();
1148  delete [] el;
1149  hcol_ = new int[bulk0_];
1150  ClpDisjointCopyN(ind,       nelems_, hcol_);
1151  mRow->nullIndexArray();
1152  delete [] ind;
1153  mrstrt_ = new CoinBigIndex[nrows_in+1];
1154  ClpDisjointCopyN(strt,  nrows_,  mrstrt_);
1155  mRow->nullStartArray();
1156  mrstrt_[nrows_] = nelems_;
1157  delete [] strt;
1158  hinrow_ = new int[nrows_in+1];
1159  ClpDisjointCopyN(len, nrows_,  hinrow_);
1160  if (nelems_>nel) {
1161    nelems_ = nel;
1162    // Clean any small elements
1163    int irow;
1164    nel=0;
1165    CoinBigIndex start=0;
1166    for (irow=0;irow<nrows_;irow++) {
1167      CoinBigIndex j;
1168      for (j=start;j<start+hinrow_[irow];j++) {
1169        hcol_[nel]=hcol_[j];
1170        if (fabs(rowels_[j])>ZTOLDP)
1171          rowels_[nel++]=rowels_[j];
1172      }
1173      start=mrstrt_[irow+1];
1174      mrstrt_[irow+1]=nel;
1175      hinrow_[irow]=nel-mrstrt_[irow];
1176    }
1177  }
1178
1179  delete mRow;
1180  if (si->integerInformation()) {
1181    CoinMemcpyN(reinterpret_cast<unsigned char *> (si->integerInformation()),ncols_,integerType_);
1182  } else {
1183    ClpFillN<unsigned char>(integerType_, ncols_, static_cast<unsigned char> (0));
1184  }
1185
1186#ifndef SLIM_CLP
1187#ifndef NO_RTTI
1188    ClpQuadraticObjective * quadraticObj = (dynamic_cast< ClpQuadraticObjective*>(si->objectiveAsObject()));
1189#else
1190    ClpQuadraticObjective * quadraticObj = NULL;
1191    if (si->objectiveAsObject()->type()==2)
1192      quadraticObj = (static_cast< ClpQuadraticObjective*>(si->objectiveAsObject()));
1193#endif
1194#endif
1195  // Set up prohibited bits if needed
1196  if (nonLinearValue) {
1197    anyProhibited_ = true;
1198    for (icol=0;icol<ncols_;icol++) {
1199      int j;
1200      bool nonLinearColumn = false;
1201      if (cost_[icol]==nonLinearValue)
1202        nonLinearColumn=true;
1203      for (j=mcstrt_[icol];j<mcstrt_[icol+1];j++) {
1204        if (colels_[j]==nonLinearValue) {
1205          nonLinearColumn=true;
1206          setRowProhibited(hrow_[j]);
1207        }
1208      }
1209      if (nonLinearColumn)
1210        setColProhibited(icol);
1211    }
1212#ifndef SLIM_CLP
1213  } else if (quadraticObj) {
1214    CoinPackedMatrix * quadratic = quadraticObj->quadraticObjective();
1215    //const int * columnQuadratic = quadratic->getIndices();
1216    //const CoinBigIndex * columnQuadraticStart = quadratic->getVectorStarts();
1217    const int * columnQuadraticLength = quadratic->getVectorLengths();
1218    //double * quadraticElement = quadratic->getMutableElements();
1219    int numberColumns = quadratic->getNumCols();
1220    anyProhibited_ = true;
1221    for (int iColumn=0;iColumn<numberColumns;iColumn++) {
1222      if (columnQuadraticLength[iColumn]) {
1223        setColProhibited(iColumn);
1224        //printf("%d prohib\n",iColumn);
1225      }
1226    }
1227#endif
1228  } else {
1229    anyProhibited_ = false;
1230  }
1231
1232  if (doStatus) {
1233    // allow for status and solution
1234    sol_ = new double[ncols_];
1235    CoinMemcpyN(si->primalColumnSolution(),ncols_,sol_);;
1236    acts_ = new double [nrows_];
1237    CoinMemcpyN(si->primalRowSolution(),nrows_,acts_);
1238    if (!si->statusArray())
1239      si->createStatus();
1240    colstat_ = new unsigned char [nrows_+ncols_];
1241    CoinMemcpyN(si->statusArray(),      (nrows_+ncols_),colstat_);
1242    rowstat_ = colstat_+ncols_;
1243  }
1244
1245  // the original model's fields are now unneeded - free them
1246 
1247  si->resize(0,0);
1248
1249#if     PRESOLVE_DEBUG
1250  matrix_bounds_ok(rlo_, rup_, nrows_);
1251  matrix_bounds_ok(clo_, cup_, ncols_);
1252#endif
1253
1254#if 0
1255  for (i=0; i<nrows; ++i)
1256    printf("NR: %6d\n", hinrow[i]);
1257  for (int i=0; i<ncols; ++i)
1258    printf("NC: %6d\n", hincol[i]);
1259#endif
1260
1261  presolve_make_memlists(mcstrt_, hincol_, clink_, ncols_);
1262  presolve_make_memlists(mrstrt_, hinrow_, rlink_, nrows_);
1263
1264  // this allows last col/row to expand up to bufsize-1 (22);
1265  // this must come after the calls to presolve_prefix
1266  mcstrt_[ncols_] = bufsize-1;
1267  mrstrt_[nrows_] = bufsize-1;
1268  // Allocate useful arrays
1269  initializeStuff();
1270
1271#if     PRESOLVE_CONSISTENCY
1272//consistent(false);
1273  presolve_consistent(this,false) ;
1274#endif
1275}
1276
1277void CoinPresolveMatrix::update_model(ClpSimplex * si,
1278                                     int nrows0,
1279                                     int ncols0,
1280                                     CoinBigIndex nelems0)
1281{
1282  si->loadProblem(ncols_, nrows_, mcstrt_, hrow_, colels_, hincol_,
1283                 clo_, cup_, cost_, rlo_, rup_);
1284  //delete [] si->integerInformation();
1285  int numberIntegers=0;
1286  for (int i=0; i<ncols_; i++) {
1287    if (integerType_[i])
1288      numberIntegers++;
1289  }
1290  if (numberIntegers) 
1291    si->copyInIntegerInformation(reinterpret_cast<const char *> (integerType_));
1292  else
1293    si->copyInIntegerInformation(NULL);
1294
1295#if     PRESOLVE_SUMMARY
1296  printf("NEW NCOL/NROW/NELS:  %d(-%d) %d(-%d) %d(-%d)\n",
1297         ncols_, ncols0-ncols_,
1298         nrows_, nrows0-nrows_,
1299         si->getNumElements(), nelems0-si->getNumElements());
1300#endif
1301  si->setDblParam(ClpObjOffset,originalOffset_-dobias_);
1302
1303}
1304
1305
1306
1307
1308
1309
1310
1311
1312
1313
1314
1315////////////////  POSTSOLVE
1316
1317CoinPostsolveMatrix::CoinPostsolveMatrix(ClpSimplex*  si,
1318                                       int ncols0_in,
1319                                       int nrows0_in,
1320                                       CoinBigIndex nelems0,
1321                                   
1322                                       double maxmin,
1323                                       // end prepost members
1324
1325                                       double *sol_in,
1326                                       double *acts_in,
1327
1328                                       unsigned char *colstat_in,
1329                                       unsigned char *rowstat_in) :
1330  CoinPrePostsolveMatrix(si,
1331                        ncols0_in, nrows0_in, nelems0,2.0),
1332
1333  free_list_(0),
1334  // link, free_list, maxlink
1335  maxlink_(bulk0_),
1336  link_(new int[/*maxlink*/ bulk0_]),
1337     
1338  cdone_(new char[ncols0_]),
1339  rdone_(new char[nrows0_in])
1340
1341{
1342  bulk0_ = maxlink_ ;
1343  nrows_ = si->getNumRows() ;
1344  ncols_ = si->getNumCols() ;
1345
1346  sol_=sol_in;
1347  rowduals_=NULL;
1348  acts_=acts_in;
1349
1350  rcosts_=NULL;
1351  colstat_=colstat_in;
1352  rowstat_=rowstat_in;
1353
1354  // this is the *reduced* model, which is probably smaller
1355  int ncols1 = ncols_ ;
1356  int nrows1 = nrows_ ;
1357
1358  const CoinPackedMatrix * m = si->matrix();
1359
1360  const CoinBigIndex nelemsr = m->getNumElements();
1361  if (m->getNumElements()&&!isGapFree(*m)) {
1362    // Odd - gaps
1363    CoinPackedMatrix mm(*m);
1364    mm.removeGaps();
1365    mm.setExtraGap(0.0);
1366   
1367    ClpDisjointCopyN(mm.getVectorStarts(), ncols1, mcstrt_);
1368    CoinZeroN(mcstrt_+ncols1,ncols0_-ncols1);
1369    mcstrt_[ncols1] = nelems0;  // ??    (should point to end of bulk store   -- lh --)
1370    ClpDisjointCopyN(mm.getVectorLengths(),ncols1,  hincol_);
1371    ClpDisjointCopyN(mm.getIndices(),      nelemsr, hrow_);
1372    ClpDisjointCopyN(mm.getElements(),     nelemsr, colels_);
1373  } else {
1374    // No gaps
1375   
1376    ClpDisjointCopyN(m->getVectorStarts(), ncols1, mcstrt_);
1377    CoinZeroN(mcstrt_+ncols1,ncols0_-ncols1);
1378    mcstrt_[ncols1] = nelems0;  // ??    (should point to end of bulk store   -- lh --)
1379    ClpDisjointCopyN(m->getVectorLengths(),ncols1,  hincol_);
1380    ClpDisjointCopyN(m->getIndices(),      nelemsr, hrow_);
1381    ClpDisjointCopyN(m->getElements(),     nelemsr, colels_);
1382  }
1383
1384
1385
1386#if     0 && PRESOLVE_DEBUG
1387  presolve_check_costs(model, &colcopy);
1388#endif
1389
1390  // This determines the size of the data structure that contains
1391  // the matrix being postsolved.  Links are taken from the free_list
1392  // to recreate matrix entries that were presolved away,
1393  // and links are added to the free_list when entries created during
1394  // presolve are discarded.  There is never a need to gc this list.
1395  // Naturally, it should contain
1396  // exactly nelems0 entries "in use" when postsolving is done,
1397  // but I don't know whether the matrix could temporarily get
1398  // larger during postsolving.  Substitution into more than two
1399  // rows could do that, in principle.  I am being very conservative
1400  // here by reserving much more than the amount of space I probably need.
1401  // If this guess is wrong, check_free_list may be called.
1402  //  int bufsize = 2*nelems0;
1403
1404  memset(cdone_, -1, ncols0_);
1405  memset(rdone_, -1, nrows0_);
1406
1407  rowduals_ = new double[nrows0_];
1408  ClpDisjointCopyN(si->getRowPrice(), nrows1, rowduals_);
1409
1410  rcosts_ = new double[ncols0_];
1411  ClpDisjointCopyN(si->getReducedCost(), ncols1, rcosts_);
1412  if (maxmin<0.0) {
1413    // change so will look as if minimize
1414    int i;
1415    for (i=0;i<nrows1;i++)
1416      rowduals_[i] = - rowduals_[i];
1417    for (i=0;i<ncols1;i++) {
1418      rcosts_[i] = - rcosts_[i];
1419    }
1420  }
1421
1422  //ClpDisjointCopyN(si->getRowUpper(), nrows1, rup_);
1423  //ClpDisjointCopyN(si->getRowLower(), nrows1, rlo_);
1424
1425  ClpDisjointCopyN(si->getColSolution(), ncols1, sol_);
1426  si->setDblParam(ClpObjOffset,originalOffset_);
1427
1428  for (int j=0; j<ncols1; j++) {
1429    CoinBigIndex kcs = mcstrt_[j];
1430    CoinBigIndex kce = kcs + hincol_[j];
1431    for (CoinBigIndex k=kcs; k<kce; ++k) {
1432      link_[k] = k+1;
1433    }
1434    link_[kce-1] = NO_LINK ;
1435  }
1436  {
1437    int ml = maxlink_;
1438    for (CoinBigIndex k=nelemsr; k<ml; ++k)
1439      link_[k] = k+1;
1440    if (ml)
1441      link_[ml-1] = NO_LINK;
1442  }
1443  free_list_ = nelemsr;
1444# if PRESOLVE_DEBUG || PRESOLVE_CONSISTENCY
1445/*
1446  These are used to track the action of postsolve transforms during debugging.
1447*/
1448  CoinFillN(cdone_,ncols1,PRESENT_IN_REDUCED) ;
1449  CoinZeroN(cdone_+ncols1,ncols0_in-ncols1) ;
1450  CoinFillN(rdone_,nrows1,PRESENT_IN_REDUCED) ;
1451  CoinZeroN(rdone_+nrows1,nrows0_in-nrows1) ;
1452# endif
1453}
1454/* This is main part of Presolve */
1455ClpSimplex * 
1456ClpPresolve::gutsOfPresolvedModel(ClpSimplex * originalModel,
1457                                  double feasibilityTolerance,
1458                                  bool keepIntegers,
1459                                  int numberPasses,
1460                                  bool dropNames,
1461                                  bool doRowObjective)
1462{
1463  ncols_ = originalModel->getNumCols();
1464  nrows_ = originalModel->getNumRows();
1465  nelems_ = originalModel->getNumElements();
1466  numberPasses_ = numberPasses;
1467
1468  double maxmin = originalModel->getObjSense();
1469  originalModel_ = originalModel;
1470  delete [] originalColumn_;
1471  originalColumn_ = new int[ncols_];
1472  delete [] originalRow_;
1473  originalRow_ = new int[nrows_];
1474  // and fill in case returns early
1475  int i;
1476  for (i=0;i<ncols_;i++)
1477    originalColumn_[i]=i;
1478  for (i=0;i<nrows_;i++)
1479    originalRow_[i]=i;
1480  delete [] rowObjective_;
1481  if (doRowObjective) {
1482    rowObjective_ = new double [nrows_];
1483    memset(rowObjective_,0,nrows_*sizeof(double));
1484  } else {
1485    rowObjective_=NULL;
1486  }
1487
1488  // result is 0 - okay, 1 infeasible, -1 go round again, 2 - original model
1489  int result = -1;
1490 
1491  // User may have deleted - its their responsibility
1492  presolvedModel_=NULL;
1493  // Messages
1494  CoinMessages messages = originalModel->coinMessages();
1495  // Only go round 100 times even if integer preprocessing
1496  int totalPasses=100;
1497  while (result==-1) {
1498
1499#ifndef CLP_NO_STD
1500    // make new copy
1501    if (saveFile_=="") {
1502#endif
1503      delete presolvedModel_;
1504#ifndef CLP_NO_STD
1505      // So won't get names
1506      int lengthNames = originalModel->lengthNames();
1507      originalModel->setLengthNames(0);
1508#endif
1509      presolvedModel_ = new ClpSimplex(*originalModel);
1510#ifndef CLP_NO_STD
1511      originalModel->setLengthNames(lengthNames);
1512    } else {
1513      presolvedModel_=originalModel;
1514    }
1515    presolvedModel_->dropNames();
1516#endif
1517
1518    // drop integer information if wanted
1519    if (!keepIntegers)
1520      presolvedModel_->deleteIntegerInformation();
1521    totalPasses--;
1522
1523    double ratio=2.0;
1524    if (substitution_>3)
1525      ratio=substitution_;
1526    else if (substitution_==2)
1527      ratio=1.5;
1528    CoinPresolveMatrix prob(ncols_,
1529                        maxmin,
1530                        presolvedModel_,
1531                        nrows_, nelems_,true,nonLinearValue_,ratio);
1532    prob.setMaximumSubstitutionLevel(substitution_);
1533    if (doRowObjective) 
1534      memset(rowObjective_,0,nrows_*sizeof(double));
1535    // See if we want statistics
1536    if ((presolveActions_&0x80000000)!=0)
1537      prob.statistics();
1538    // make sure row solution correct
1539    {
1540      double *colels    = prob.colels_;
1541      int *hrow         = prob.hrow_;
1542      CoinBigIndex *mcstrt              = prob.mcstrt_;
1543      int *hincol               = prob.hincol_;
1544      int ncols         = prob.ncols_;
1545     
1546     
1547      double * csol = prob.sol_;
1548      double * acts = prob.acts_;
1549      int nrows = prob.nrows_;
1550
1551      int colx;
1552
1553      memset(acts,0,nrows*sizeof(double));
1554     
1555      for (colx = 0; colx < ncols; ++colx) {
1556        double solutionValue = csol[colx];
1557        for (int i=mcstrt[colx]; i<mcstrt[colx]+hincol[colx]; ++i) {
1558          int row = hrow[i];
1559          double coeff = colels[i];
1560          acts[row] += solutionValue*coeff;
1561        }
1562      }
1563    }
1564
1565    // move across feasibility tolerance
1566    prob.feasibilityTolerance_ = feasibilityTolerance;
1567
1568    // Do presolve
1569    paction_ = presolve(&prob);
1570    // Get rid of useful arrays
1571    prob.deleteStuff();
1572
1573    result =0; 
1574
1575    if (prob.status_==0&&paction_) {
1576      // Looks feasible but double check to see if anything slipped through
1577      int n             = prob.ncols_;
1578      double * lo = prob.clo_;
1579      double * up = prob.cup_;
1580      int i;
1581     
1582      for (i=0;i<n;i++) {
1583        if (up[i]<lo[i]) {
1584          if (up[i]<lo[i]-1.0e-8) {
1585            // infeasible
1586            prob.status_=1;
1587          } else {
1588            up[i]=lo[i];
1589          }
1590        }
1591      }
1592     
1593      n = prob.nrows_;
1594      lo = prob.rlo_;
1595      up = prob.rup_;
1596
1597      for (i=0;i<n;i++) {
1598        if (up[i]<lo[i]) {
1599          if (up[i]<lo[i]-1.0e-8) {
1600            // infeasible
1601            prob.status_=1;
1602          } else {
1603            up[i]=lo[i];
1604          }
1605        }
1606      }
1607    }
1608    if (prob.status_==0&&paction_) {
1609      // feasible
1610   
1611      prob.update_model(presolvedModel_, nrows_, ncols_, nelems_);
1612      // copy status and solution
1613      CoinMemcpyN(           prob.sol_,prob.ncols_,presolvedModel_->primalColumnSolution());
1614      CoinMemcpyN(           prob.acts_,prob.nrows_,presolvedModel_->primalRowSolution());
1615      CoinMemcpyN(           prob.colstat_,prob.ncols_,presolvedModel_->statusArray());
1616      CoinMemcpyN(           prob.rowstat_,prob.nrows_,presolvedModel_->statusArray()+prob.ncols_);
1617      delete [] prob.sol_;
1618      delete [] prob.acts_;
1619      delete [] prob.colstat_;
1620      prob.sol_=NULL;
1621      prob.acts_=NULL;
1622      prob.colstat_=NULL;
1623     
1624      int ncolsNow = presolvedModel_->getNumCols();
1625      CoinMemcpyN(prob.originalColumn_,ncolsNow,originalColumn_);
1626#ifndef SLIM_CLP
1627#ifndef NO_RTTI
1628      ClpQuadraticObjective * quadraticObj = (dynamic_cast< ClpQuadraticObjective*>(originalModel->objectiveAsObject()));
1629#else
1630      ClpQuadraticObjective * quadraticObj = NULL;
1631      if (originalModel->objectiveAsObject()->type()==2)
1632        quadraticObj = (static_cast< ClpQuadraticObjective*>(originalModel->objectiveAsObject()));
1633#endif
1634      if (quadraticObj) {
1635        // set up for subset
1636        char * mark = new char [ncols_];
1637        memset(mark,0,ncols_);
1638        CoinPackedMatrix * quadratic = quadraticObj->quadraticObjective();
1639        //const int * columnQuadratic = quadratic->getIndices();
1640        //const CoinBigIndex * columnQuadraticStart = quadratic->getVectorStarts();
1641        const int * columnQuadraticLength = quadratic->getVectorLengths();
1642        //double * quadraticElement = quadratic->getMutableElements();
1643        int numberColumns = quadratic->getNumCols();
1644        ClpQuadraticObjective * newObj = new ClpQuadraticObjective(*quadraticObj,
1645                                                                   ncolsNow,
1646                                                                   originalColumn_);
1647        // and modify linear and check
1648        double * linear = newObj->linearObjective();
1649 CoinMemcpyN(presolvedModel_->objective(),ncolsNow,linear);
1650        int iColumn;
1651        for ( iColumn=0;iColumn<numberColumns;iColumn++) {
1652          if (columnQuadraticLength[iColumn])
1653            mark[iColumn]=1;
1654        }
1655        // and new
1656        quadratic = newObj->quadraticObjective();
1657        columnQuadraticLength = quadratic->getVectorLengths();
1658        int numberColumns2 = quadratic->getNumCols();
1659        for ( iColumn=0;iColumn<numberColumns2;iColumn++) {
1660          if (columnQuadraticLength[iColumn])
1661            mark[originalColumn_[iColumn]]=0;
1662        }
1663        presolvedModel_->setObjective(newObj);
1664        delete newObj;
1665        // final check
1666        for ( iColumn=0;iColumn<numberColumns;iColumn++) 
1667          if (mark[iColumn])
1668            printf("Quadratic column %d modified - may be okay\n",iColumn);
1669        delete [] mark;
1670      }
1671#endif
1672      delete [] prob.originalColumn_;
1673      prob.originalColumn_=NULL;
1674      int nrowsNow = presolvedModel_->getNumRows();
1675      CoinMemcpyN(prob.originalRow_,nrowsNow,originalRow_);
1676      delete [] prob.originalRow_;
1677      prob.originalRow_=NULL;
1678#ifndef CLP_NO_STD
1679      if (!dropNames&&originalModel->lengthNames()) {
1680        // Redo names
1681        int iRow;
1682        std::vector<std::string> rowNames;
1683        rowNames.reserve(nrowsNow);
1684        for (iRow=0;iRow<nrowsNow;iRow++) {
1685          int kRow = originalRow_[iRow];
1686          rowNames.push_back(originalModel->rowName(kRow));
1687        }
1688     
1689        int iColumn;
1690        std::vector<std::string> columnNames;
1691        columnNames.reserve(ncolsNow);
1692        for (iColumn=0;iColumn<ncolsNow;iColumn++) {
1693          int kColumn = originalColumn_[iColumn];
1694          columnNames.push_back(originalModel->columnName(kColumn));
1695        }
1696        presolvedModel_->copyNames(rowNames,columnNames);
1697      } else {
1698        presolvedModel_->setLengthNames(0);
1699      }
1700#endif
1701      if (rowObjective_) {
1702        int iRow;
1703        int k=-1;
1704        int nObj=0;
1705        for (iRow=0;iRow<nrowsNow;iRow++) {
1706          int kRow = originalRow_[iRow];
1707          assert (kRow>k);
1708          k=kRow;
1709          rowObjective_[iRow]=rowObjective_[kRow];
1710          if (rowObjective_[iRow])
1711            nObj++;
1712        }
1713        if (nObj) {
1714          printf("%d costed slacks\n",nObj);
1715          presolvedModel_->setRowObjective(rowObjective_);
1716        }
1717      }
1718      // now clean up integer variables.  This can modify original
1719      int i;
1720      const char * information = presolvedModel_->integerInformation();
1721      if (information) {
1722        int numberChanges=0;
1723        double * lower0 = originalModel_->columnLower();
1724        double * upper0 = originalModel_->columnUpper();
1725        double * lower = presolvedModel_->columnLower();
1726        double * upper = presolvedModel_->columnUpper();
1727        for (i=0;i<ncolsNow;i++) {
1728          if (!information[i])
1729            continue;
1730          int iOriginal = originalColumn_[i];
1731          double lowerValue0 = lower0[iOriginal];
1732          double upperValue0 = upper0[iOriginal];
1733          double lowerValue = ceil(lower[i]-1.0e-5);
1734          double upperValue = floor(upper[i]+1.0e-5);
1735          lower[i]=lowerValue;
1736          upper[i]=upperValue;
1737          if (lowerValue>upperValue) {
1738            numberChanges++;
1739            presolvedModel_->messageHandler()->message(COIN_PRESOLVE_COLINFEAS,
1740                                                       messages)
1741                                                         <<iOriginal
1742                                                         <<lowerValue
1743                                                         <<upperValue
1744                                                         <<CoinMessageEol;
1745            result=1;
1746          } else {
1747            if (lowerValue>lowerValue0+1.0e-8) {
1748              lower0[iOriginal] = lowerValue;
1749              numberChanges++;
1750            }
1751            if (upperValue<upperValue0-1.0e-8) {
1752              upper0[iOriginal] = upperValue;
1753              numberChanges++;
1754            }
1755          }       
1756        }
1757        if (numberChanges) {
1758          presolvedModel_->messageHandler()->message(COIN_PRESOLVE_INTEGERMODS,
1759                                                     messages)
1760                                                       <<numberChanges
1761                                                       <<CoinMessageEol;
1762          if (!result&&totalPasses>0) {
1763            result = -1; // round again
1764            const CoinPresolveAction *paction = paction_;
1765            while (paction) {
1766              const CoinPresolveAction *next = paction->next;
1767              delete paction;
1768              paction = next;
1769            }
1770            paction_=NULL;
1771          }
1772        }
1773      }
1774    } else if (prob.status_) {
1775      // infeasible or unbounded
1776      result=1;
1777      originalModel->setProblemStatus(prob.status_);
1778    } else {
1779      // no changes - model needs restoring after Lou's changes
1780#ifndef CLP_NO_STD
1781      if (saveFile_=="") {
1782#endif
1783        delete presolvedModel_;
1784        presolvedModel_ = new ClpSimplex(*originalModel);
1785        // but we need to remove gaps
1786        ClpPackedMatrix* clpMatrix =
1787          dynamic_cast< ClpPackedMatrix*>(presolvedModel_->clpMatrix());
1788        if (clpMatrix) {
1789          clpMatrix->getPackedMatrix()->removeGaps();
1790        }
1791#ifndef CLP_NO_STD
1792      } else {
1793        presolvedModel_=originalModel;
1794      }
1795      presolvedModel_->dropNames();
1796#endif
1797     
1798      // drop integer information if wanted
1799      if (!keepIntegers)
1800        presolvedModel_->deleteIntegerInformation();
1801      result=2;
1802    }
1803  }
1804  if (result==0||result==2) {
1805    int nrowsAfter = presolvedModel_->getNumRows();
1806    int ncolsAfter = presolvedModel_->getNumCols();
1807    CoinBigIndex nelsAfter = presolvedModel_->getNumElements();
1808    presolvedModel_->messageHandler()->message(COIN_PRESOLVE_STATS,
1809                                               messages)
1810                                                 <<nrowsAfter<< -(nrows_ - nrowsAfter)
1811                                                 << ncolsAfter<< -(ncols_ - ncolsAfter)
1812                                                 <<nelsAfter<< -(nelems_ - nelsAfter)
1813                                                 <<CoinMessageEol;
1814  } else {
1815    destroyPresolve();
1816    if (presolvedModel_!=originalModel_)
1817      delete presolvedModel_;
1818    presolvedModel_=NULL;
1819  }
1820  return presolvedModel_;
1821}
1822
1823
Note: See TracBrowser for help on using the repository browser.