source: branches/devel/Clp/src/ClpPresolve.cpp @ 867

Last change on this file since 867 was 867, checked in by forrest, 14 years ago

integers in save/restore
allow options in Coinsolve

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