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

Last change on this file since 1051 was 1051, checked in by forrest, 12 years ago

to eliminate memory leak

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