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

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

many changes to try and improve performance

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 47.8 KB
Line 
1// Copyright (C) 2002, International Business Machines
2// Corporation and others.  All Rights Reserved.
3
4//#define       PRESOLVE_CONSISTENCY    1
5//#define       PRESOLVE_DEBUG  1
6
7#include <stdio.h>
8
9#include <assert.h>
10#include <iostream>
11
12#include "CoinHelperFunctions.hpp"
13
14#include "CoinPackedMatrix.hpp"
15#include "ClpPackedMatrix.hpp"
16#include "ClpSimplex.hpp"
17#ifndef SLIM_CLP
18#include "ClpQuadraticObjective.hpp"
19#endif
20
21#include "ClpPresolve.hpp"
22#include "CoinPresolveMatrix.hpp"
23
24#include "CoinPresolveEmpty.hpp"
25#include "CoinPresolveFixed.hpp"
26#include "CoinPresolvePsdebug.hpp"
27#include "CoinPresolveSingleton.hpp"
28#include "CoinPresolveDoubleton.hpp"
29#include "CoinPresolveTripleton.hpp"
30#include "CoinPresolveZeros.hpp"
31#include "CoinPresolveSubst.hpp"
32#include "CoinPresolveForcing.hpp"
33#include "CoinPresolveDual.hpp"
34#include "CoinPresolveTighten.hpp"
35#include "CoinPresolveUseless.hpp"
36#include "CoinPresolveDupcol.hpp"
37#include "CoinPresolveImpliedFree.hpp"
38#include "CoinPresolveIsolated.hpp"
39#include "CoinMessage.hpp"
40
41
42
43ClpPresolve::ClpPresolve() :
44  originalModel_(NULL),
45  presolvedModel_(NULL),
46  nonLinearValue_(0.0),
47  originalColumn_(NULL),
48  originalRow_(NULL),
49  rowObjective_(NULL),
50  paction_(0),   
51  ncols_(0),
52  nrows_(0),
53  nelems_(0),
54  numberPasses_(5),
55  substitution_(3),
56#ifndef CLP_NO_STD
57  saveFile_(""),
58#endif
59  presolveActions_(0)
60{
61}
62
63ClpPresolve::~ClpPresolve()
64{
65  destroyPresolve();
66}
67// Gets rid of presolve actions (e.g.when infeasible)
68void 
69ClpPresolve::destroyPresolve()
70{
71  const CoinPresolveAction *paction = paction_;
72  while (paction) {
73    const CoinPresolveAction *next = paction->next;
74    delete paction;
75    paction = next;
76  }
77  delete [] originalColumn_;
78  delete [] originalRow_;
79  paction_=NULL;
80  originalColumn_=NULL;
81  originalRow_=NULL;
82  delete [] rowObjective_;
83  rowObjective_=NULL;
84}
85
86/* This version of presolve returns a pointer to a new presolved
87   model.  NULL if infeasible
88*/
89ClpSimplex * 
90ClpPresolve::presolvedModel(ClpSimplex & si,
91                            double feasibilityTolerance,
92                            bool keepIntegers,
93                            int numberPasses,
94                            bool dropNames,
95                            bool doRowObjective)
96{
97  // Check matrix
98  if (!si.clpMatrix()->allElementsInRange(&si,si.getSmallElementValue(),
99                                          1.0e20))
100    return NULL;
101  else
102    return gutsOfPresolvedModel(&si,feasibilityTolerance,keepIntegers,numberPasses,dropNames,
103                                doRowObjective);
104}
105#ifndef CLP_NO_STD
106/* This version of presolve updates
107   model and saves original data to file.  Returns non-zero if infeasible
108*/
109int
110ClpPresolve::presolvedModelToFile(ClpSimplex &si,std::string fileName,
111                                  double feasibilityTolerance,
112                                  bool keepIntegers,
113                                  int numberPasses,
114                                  bool doRowObjective)
115{
116  // Check matrix
117  if (!si.clpMatrix()->allElementsInRange(&si,si.getSmallElementValue(),
118                                          1.0e20))
119    return 2;
120  saveFile_=fileName;
121  si.saveModel(saveFile_.c_str());
122  ClpSimplex * model = gutsOfPresolvedModel(&si,feasibilityTolerance,keepIntegers,numberPasses,true,
123                                            doRowObjective);
124  if (model==&si) {
125    return 0;
126  } else {
127    si.restoreModel(saveFile_.c_str());
128    remove(saveFile_.c_str());
129    return 1;
130  }
131}
132#endif
133// Return pointer to presolved model
134ClpSimplex * 
135ClpPresolve::model() const
136{
137  return presolvedModel_;
138}
139// Return pointer to original model
140ClpSimplex * 
141ClpPresolve::originalModel() const
142{
143  return originalModel_;
144}
145void 
146ClpPresolve::postsolve(bool updateStatus)
147{
148  // Return at once if no presolved model
149  if (!presolvedModel_)
150    return;
151  // Messages
152  CoinMessages messages = originalModel_->coinMessages();
153  if (!presolvedModel_->isProvenOptimal()) {
154    presolvedModel_->messageHandler()->message(COIN_PRESOLVE_NONOPTIMAL,
155                                             messages)
156                                               <<CoinMessageEol;
157  }
158
159  // this is the size of the original problem
160  const int ncols0  = ncols_;
161  const int nrows0  = nrows_;
162  const CoinBigIndex nelems0 = nelems_;
163 
164  // this is the reduced problem
165  int ncols = presolvedModel_->getNumCols();
166  int nrows = presolvedModel_->getNumRows();
167
168  double * acts=NULL;
169  double * sol =NULL;
170  unsigned char * rowstat=NULL;
171  unsigned char * colstat = NULL;
172#ifndef CLP_NO_STD
173  if (saveFile_=="") {
174#endif
175    // reality check
176    assert(ncols0==originalModel_->getNumCols());
177    assert(nrows0==originalModel_->getNumRows());
178    acts = originalModel_->primalRowSolution();
179    sol  = originalModel_->primalColumnSolution();
180    if (updateStatus) {
181      // postsolve does not know about fixed
182      int i;
183      for (i=0;i<nrows+ncols;i++) {
184        if (presolvedModel_->getColumnStatus(i)==ClpSimplex::isFixed)
185          presolvedModel_->setColumnStatus(i,ClpSimplex::atLowerBound);
186      }
187      unsigned char *status = originalModel_->statusArray();
188      if (!status) {
189        originalModel_->createStatus();
190        status = originalModel_->statusArray();
191      }
192      rowstat = status + ncols0;
193      colstat = status;
194      CoinMemcpyN( presolvedModel_->statusArray(), ncols,colstat);
195      CoinMemcpyN( presolvedModel_->statusArray()+ncols, nrows,rowstat);
196    }
197#ifndef CLP_NO_STD
198  } else {
199    // from file
200    acts = new double[nrows0];
201    sol  = new double[ncols0];
202    CoinZeroN(acts,nrows0);
203    CoinZeroN(sol,ncols0);
204    if (updateStatus) {
205      unsigned char *status = new unsigned char [nrows0+ncols0];
206      rowstat = status + ncols0;
207      colstat = status;
208      CoinMemcpyN( presolvedModel_->statusArray(), ncols,colstat);
209      CoinMemcpyN( presolvedModel_->statusArray()+ncols, nrows,rowstat);
210    }
211  }
212#endif
213
214  // CoinPostsolveMatrix object assumes ownership of sol, acts, colstat;
215  // will be deleted by ~CoinPostsolveMatrix. delete[] operations below
216  // cause duplicate free. In case where saveFile == "", as best I can see
217  // arrays are owned by originalModel_. fix is to
218  // clear fields in prob to avoid delete[] in ~CoinPostsolveMatrix.
219  CoinPostsolveMatrix prob(presolvedModel_,
220                       ncols0,
221                       nrows0,
222                       nelems0,
223                       presolvedModel_->getObjSense(),
224                       // end prepost
225                       
226                       sol, acts,
227                       colstat, rowstat);
228   
229  postsolve(prob);
230
231#ifndef CLP_NO_STD
232  if (saveFile_!="") {
233    // From file
234    assert (originalModel_==presolvedModel_);
235    originalModel_->restoreModel(saveFile_.c_str());
236    remove(saveFile_.c_str());
237    CoinMemcpyN(acts,nrows0,originalModel_->primalRowSolution());
238    // delete [] acts;
239    CoinMemcpyN(sol,ncols0,originalModel_->primalColumnSolution());
240    // delete [] sol;
241    if (updateStatus) {
242      CoinMemcpyN(colstat,nrows0+ncols0,originalModel_->statusArray());
243      // delete [] colstat;
244    }
245  } else {
246#endif
247    prob.sol_ = 0 ;
248    prob.acts_ = 0 ;
249    prob.colstat_ = 0 ;
250#ifndef CLP_NO_STD
251  }
252#endif
253  // put back duals
254  CoinMemcpyN(prob.rowduals_,   nrows_,originalModel_->dualRowSolution());
255  double maxmin = originalModel_->getObjSense();
256  if (maxmin<0.0) {
257    // swap signs
258    int i;
259    double * pi = originalModel_->dualRowSolution();
260    for (i=0;i<nrows_;i++)
261      pi[i] = -pi[i];
262  }
263  // Now check solution
264  double offset;
265  CoinMemcpyN(originalModel_->objectiveAsObject()->gradient(originalModel_,
266                                                            originalModel_->primalColumnSolution(),offset,true),
267              ncols_,originalModel_->dualColumnSolution());
268  originalModel_->transposeTimes(-1.0,
269                                 originalModel_->dualRowSolution(),
270                                 originalModel_->dualColumnSolution());
271  memset(originalModel_->primalRowSolution(),0,nrows_*sizeof(double));
272  originalModel_->times(1.0,originalModel_->primalColumnSolution(),
273                        originalModel_->primalRowSolution());
274  originalModel_->checkSolutionInternal();
275  // Messages
276  presolvedModel_->messageHandler()->message(COIN_PRESOLVE_POSTSOLVE,
277                                            messages)
278                                              <<originalModel_->objectiveValue()
279                                              <<originalModel_->sumDualInfeasibilities()
280                                              <<originalModel_->numberDualInfeasibilities()
281                                              <<originalModel_->sumPrimalInfeasibilities()
282                                              <<originalModel_->numberPrimalInfeasibilities()
283                                               <<CoinMessageEol;
284 
285  //originalModel_->objectiveValue_=objectiveValue_;
286  originalModel_->setNumberIterations(presolvedModel_->numberIterations());
287  if (!presolvedModel_->status()) {
288    if (!originalModel_->numberDualInfeasibilities()&&
289        !originalModel_->numberPrimalInfeasibilities()) {
290      originalModel_->setProblemStatus( 0);
291    } else {
292      originalModel_->setProblemStatus( -1);
293      // Say not optimal after presolve
294      originalModel_->setSecondaryStatus(7);
295      presolvedModel_->messageHandler()->message(COIN_PRESOLVE_NEEDS_CLEANING,
296                                            messages)
297                                              <<CoinMessageEol;
298    }
299  } else {
300    originalModel_->setProblemStatus( presolvedModel_->status());
301  }
302#ifndef CLP_NO_STD
303  if (saveFile_!="") 
304    presolvedModel_=NULL;
305#endif
306}
307
308// return pointer to original columns
309const int * 
310ClpPresolve::originalColumns() const
311{
312  return originalColumn_;
313}
314// return pointer to original rows
315const int * 
316ClpPresolve::originalRows() const
317{
318  return originalRow_;
319}
320// Set pointer to original model
321void 
322ClpPresolve::setOriginalModel(ClpSimplex * model)
323{
324  originalModel_=model;
325}
326#if 0
327// A lazy way to restrict which transformations are applied
328// during debugging.
329static int ATOI(const char *name)
330{
331 return true;
332#if     PRESOLVE_DEBUG || PRESOLVE_SUMMARY
333  if (getenv(name)) {
334    int val = atoi(getenv(name));
335    printf("%s = %d\n", name, val);
336    return (val);
337  } else {
338    if (strcmp(name,"off"))
339      return (true);
340    else
341      return (false);
342  }
343#else
344  return (true);
345#endif
346}
347#endif
348//#define PRESOLVE_DEBUG 1
349#if PRESOLVE_DEBUG
350void check_sol(CoinPresolveMatrix *prob,double tol)
351{
352  double *colels        = prob->colels_;
353  int *hrow             = prob->hrow_;
354  int *mcstrt           = prob->mcstrt_;
355  int *hincol           = prob->hincol_;
356  int *hinrow           = prob->hinrow_;
357  int ncols             = prob->ncols_;
358
359
360  double * csol = prob->sol_;
361  double * acts = prob->acts_;
362  double * clo = prob->clo_;
363  double * cup = prob->cup_;
364  int nrows = prob->nrows_;
365  double * rlo = prob->rlo_;
366  double * rup = prob->rup_;
367
368  int colx;
369
370  double * rsol = new double[nrows];
371  memset(rsol,0,nrows*sizeof(double));
372
373  for (colx = 0; colx < ncols; ++colx) {
374    if (1) {
375      CoinBigIndex k = mcstrt[colx];
376      int nx = hincol[colx];
377      double solutionValue = csol[colx];
378      for (int i=0; i<nx; ++i) {
379        int row = hrow[k];
380        double coeff = colels[k];
381        k++;
382        rsol[row] += solutionValue*coeff;
383      }
384      if (csol[colx]<clo[colx]-tol) {
385        printf("low CSOL:  %d  - %g %g %g\n",
386                   colx, clo[colx], csol[colx], cup[colx]);
387      } else if (csol[colx]>cup[colx]+tol) {
388        printf("high CSOL:  %d  - %g %g %g\n",
389                   colx, clo[colx], csol[colx], cup[colx]);
390      } 
391    }
392  }
393  int rowx;
394  for (rowx = 0; rowx < nrows; ++rowx) {
395    if (hinrow[rowx]) {
396      if (fabs(rsol[rowx]-acts[rowx])>tol)
397        printf("inacc RSOL:  %d - %g %g (acts_ %g) %g\n",
398                   rowx,  rlo[rowx], rsol[rowx], acts[rowx], rup[rowx]);
399      if (rsol[rowx]<rlo[rowx]-tol) {
400        printf("low RSOL:  %d - %g %g %g\n",
401                   rowx,  rlo[rowx], rsol[rowx], rup[rowx]);
402      } else if (rsol[rowx]>rup[rowx]+tol ) {
403        printf("high RSOL:  %d - %g %g %g\n",
404                   rowx,  rlo[rowx], rsol[rowx], rup[rowx]);
405      } 
406    }
407  }
408  delete [] rsol;
409}
410#endif
411// This is the presolve loop.
412// It is a separate virtual function so that it can be easily
413// customized by subclassing CoinPresolve.
414const CoinPresolveAction *ClpPresolve::presolve(CoinPresolveMatrix *prob)
415{
416  // Messages
417  CoinMessages messages = CoinMessage(prob->messages().language());
418  paction_ = 0;
419
420  prob->status_=0; // say feasible
421  paction_ = make_fixed(prob, paction_);
422  // if integers then switch off dual stuff
423  // later just do individually
424  bool doDualStuff = (presolvedModel_->integerInformation()==NULL);
425  // but allow in some cases
426  if ((presolveActions_&512)!=0)
427    doDualStuff=true;
428  if (prob->anyProhibited()) 
429    doDualStuff=false;
430  if (!doDual())
431    doDualStuff=false;
432#if     PRESOLVE_CONSISTENCY
433//  presolve_links_ok(prob->rlink_, prob->mrstrt_, prob->hinrow_, prob->nrows_);
434    presolve_links_ok(prob,false,true) ;
435#endif
436
437  if (!prob->status_) {
438    bool slackSingleton = doSingletonColumn();
439    slackSingleton = true;
440    const bool slackd = doSingleton();
441    const bool doubleton = doDoubleton();
442    const bool tripleton = doTripleton();
443    const bool forcing = doForcing();
444    const bool ifree = doImpliedFree();
445    const bool zerocost = doTighten();
446    const bool dupcol = doDupcol();
447    const bool duprow = doDuprow();
448    const bool dual = doDualStuff;
449   
450    // some things are expensive so just do once (normally)
451
452    int i;
453    // say look at all
454    if (!prob->anyProhibited()) {
455      for (i=0;i<nrows_;i++) 
456        prob->rowsToDo_[i]=i;
457      prob->numberRowsToDo_=nrows_;
458      for (i=0;i<ncols_;i++) 
459        prob->colsToDo_[i]=i;
460      prob->numberColsToDo_=ncols_;
461    } else {
462      // some stuff must be left alone
463      prob->numberRowsToDo_=0;
464      for (i=0;i<nrows_;i++) 
465        if (!prob->rowProhibited(i))
466            prob->rowsToDo_[prob->numberRowsToDo_++]=i;
467      prob->numberColsToDo_=0;
468      for (i=0;i<ncols_;i++) 
469        if (!prob->colProhibited(i))
470            prob->colsToDo_[prob->numberColsToDo_++]=i;
471    }
472
473
474    int iLoop;
475#if     PRESOLVE_DEBUG
476    check_sol(prob,1.0e0);
477#endif
478    if (dupcol) {
479      //paction_ = dupcol_action::presolve(prob, paction_);
480    }
481
482    if ((presolveActions_&16384)!=0)
483      prob->setPresolveOptions(prob->presolveOptions()|16384);
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    CoinMemcpyN((unsigned char *) si->integerInformation(),ncols_,integerType_);
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    CoinMemcpyN(si->primalColumnSolution(),ncols_,sol_);;
1177    acts_ = new double [nrows_];
1178    CoinMemcpyN(si->primalRowSolution(),nrows_,acts_);
1179    if (!si->statusArray())
1180      si->createStatus();
1181    colstat_ = new unsigned char [nrows_+ncols_];
1182    CoinMemcpyN(si->statusArray(),      (nrows_+ncols_),colstat_);
1183    rowstat_ = colstat_+ncols_;
1184  }
1185
1186  // the original model's fields are now unneeded - free them
1187 
1188  si->resize(0,0);
1189
1190#if     PRESOLVE_DEBUG
1191  matrix_bounds_ok(rlo_, rup_, nrows_);
1192  matrix_bounds_ok(clo_, cup_, ncols_);
1193#endif
1194
1195#if 0
1196  for (i=0; i<nrows; ++i)
1197    printf("NR: %6d\n", hinrow[i]);
1198  for (int i=0; i<ncols; ++i)
1199    printf("NC: %6d\n", hincol[i]);
1200#endif
1201
1202  presolve_make_memlists(mcstrt_, hincol_, clink_, ncols_);
1203  presolve_make_memlists(mrstrt_, hinrow_, rlink_, nrows_);
1204
1205  // this allows last col/row to expand up to bufsize-1 (22);
1206  // this must come after the calls to presolve_prefix
1207  mcstrt_[ncols_] = bufsize-1;
1208  mrstrt_[nrows_] = bufsize-1;
1209
1210#if     PRESOLVE_CONSISTENCY
1211//consistent(false);
1212  presolve_consistent(this,false) ;
1213#endif
1214}
1215
1216void CoinPresolveMatrix::update_model(ClpSimplex * si,
1217                                     int nrows0,
1218                                     int ncols0,
1219                                     CoinBigIndex nelems0)
1220{
1221  si->loadProblem(ncols_, nrows_, mcstrt_, hrow_, colels_, hincol_,
1222                 clo_, cup_, cost_, rlo_, rup_);
1223  //delete [] si->integerInformation();
1224  int numberIntegers=0;
1225  for (int i=0; i<ncols_; i++) {
1226    if (integerType_[i])
1227      numberIntegers++;
1228  }
1229  if (numberIntegers) 
1230    si->copyInIntegerInformation((const char *) integerType_);
1231  else
1232    si->copyInIntegerInformation(NULL);
1233
1234#if     PRESOLVE_SUMMARY
1235  printf("NEW NCOL/NROW/NELS:  %d(-%d) %d(-%d) %d(-%d)\n",
1236         ncols_, ncols0-ncols_,
1237         nrows_, nrows0-nrows_,
1238         si->getNumElements(), nelems0-si->getNumElements());
1239#endif
1240  si->setDblParam(ClpObjOffset,originalOffset_-dobias_);
1241
1242}
1243
1244
1245
1246
1247
1248
1249
1250
1251
1252
1253
1254////////////////  POSTSOLVE
1255
1256CoinPostsolveMatrix::CoinPostsolveMatrix(ClpSimplex*  si,
1257                                       int ncols0_in,
1258                                       int nrows0_in,
1259                                       CoinBigIndex nelems0,
1260                                   
1261                                       double maxmin,
1262                                       // end prepost members
1263
1264                                       double *sol_in,
1265                                       double *acts_in,
1266
1267                                       unsigned char *colstat_in,
1268                                       unsigned char *rowstat_in) :
1269  CoinPrePostsolveMatrix(si,
1270                        ncols0_in, nrows0_in, nelems0,2.0),
1271
1272  free_list_(0),
1273  // link, free_list, maxlink
1274  maxlink_(bulk0_),
1275  link_(new int[/*maxlink*/ bulk0_]),
1276     
1277  cdone_(new char[ncols0_]),
1278  rdone_(new char[nrows0_in])
1279
1280{
1281  bulk0_ = maxlink_ ;
1282  nrows_ = si->getNumRows() ;
1283  ncols_ = si->getNumCols() ;
1284
1285  sol_=sol_in;
1286  rowduals_=NULL;
1287  acts_=acts_in;
1288
1289  rcosts_=NULL;
1290  colstat_=colstat_in;
1291  rowstat_=rowstat_in;
1292
1293  // this is the *reduced* model, which is probably smaller
1294  int ncols1 = ncols_ ;
1295  int nrows1 = nrows_ ;
1296
1297  const CoinPackedMatrix * m = si->matrix();
1298
1299  const CoinBigIndex nelemsr = m->getNumElements();
1300  if (m->getNumElements()&&!isGapFree(*m)) {
1301    // Odd - gaps
1302    CoinPackedMatrix mm(*m);
1303    mm.removeGaps();
1304    mm.setExtraGap(0.0);
1305   
1306    ClpDisjointCopyN(mm.getVectorStarts(), ncols1, mcstrt_);
1307    CoinZeroN(mcstrt_+ncols1,ncols0_-ncols1);
1308    mcstrt_[ncols1] = nelems0;  // ??    (should point to end of bulk store   -- lh --)
1309    ClpDisjointCopyN(mm.getVectorLengths(),ncols1,  hincol_);
1310    ClpDisjointCopyN(mm.getIndices(),      nelemsr, hrow_);
1311    ClpDisjointCopyN(mm.getElements(),     nelemsr, colels_);
1312  } else {
1313    // No gaps
1314   
1315    ClpDisjointCopyN(m->getVectorStarts(), ncols1, mcstrt_);
1316    CoinZeroN(mcstrt_+ncols1,ncols0_-ncols1);
1317    mcstrt_[ncols1] = nelems0;  // ??    (should point to end of bulk store   -- lh --)
1318    ClpDisjointCopyN(m->getVectorLengths(),ncols1,  hincol_);
1319    ClpDisjointCopyN(m->getIndices(),      nelemsr, hrow_);
1320    ClpDisjointCopyN(m->getElements(),     nelemsr, colels_);
1321  }
1322
1323
1324
1325#if     0 && PRESOLVE_DEBUG
1326  presolve_check_costs(model, &colcopy);
1327#endif
1328
1329  // This determines the size of the data structure that contains
1330  // the matrix being postsolved.  Links are taken from the free_list
1331  // to recreate matrix entries that were presolved away,
1332  // and links are added to the free_list when entries created during
1333  // presolve are discarded.  There is never a need to gc this list.
1334  // Naturally, it should contain
1335  // exactly nelems0 entries "in use" when postsolving is done,
1336  // but I don't know whether the matrix could temporarily get
1337  // larger during postsolving.  Substitution into more than two
1338  // rows could do that, in principle.  I am being very conservative
1339  // here by reserving much more than the amount of space I probably need.
1340  // If this guess is wrong, check_free_list may be called.
1341  //  int bufsize = 2*nelems0;
1342
1343  memset(cdone_, -1, ncols0_);
1344  memset(rdone_, -1, nrows0_);
1345
1346  rowduals_ = new double[nrows0_];
1347  ClpDisjointCopyN(si->getRowPrice(), nrows1, rowduals_);
1348
1349  rcosts_ = new double[ncols0_];
1350  ClpDisjointCopyN(si->getReducedCost(), ncols1, rcosts_);
1351  if (maxmin<0.0) {
1352    // change so will look as if minimize
1353    int i;
1354    for (i=0;i<nrows1;i++)
1355      rowduals_[i] = - rowduals_[i];
1356    for (i=0;i<ncols1;i++) {
1357      rcosts_[i] = - rcosts_[i];
1358    }
1359  }
1360
1361  //ClpDisjointCopyN(si->getRowUpper(), nrows1, rup_);
1362  //ClpDisjointCopyN(si->getRowLower(), nrows1, rlo_);
1363
1364  ClpDisjointCopyN(si->getColSolution(), ncols1, sol_);
1365  si->setDblParam(ClpObjOffset,originalOffset_);
1366
1367  for (int j=0; j<ncols1; j++) {
1368    CoinBigIndex kcs = mcstrt_[j];
1369    CoinBigIndex kce = kcs + hincol_[j];
1370    for (CoinBigIndex k=kcs; k<kce; ++k) {
1371      link_[k] = k+1;
1372    }
1373    link_[kce-1] = NO_LINK ;
1374  }
1375  {
1376    int ml = maxlink_;
1377    for (CoinBigIndex k=nelemsr; k<ml; ++k)
1378      link_[k] = k+1;
1379    if (ml)
1380      link_[ml-1] = NO_LINK;
1381  }
1382  free_list_ = nelemsr;
1383# if PRESOLVE_DEBUG || PRESOLVE_CONSISTENCY
1384/*
1385  These are used to track the action of postsolve transforms during debugging.
1386*/
1387  CoinFillN(cdone_,ncols1,PRESENT_IN_REDUCED) ;
1388  CoinZeroN(cdone_+ncols1,ncols0_in-ncols1) ;
1389  CoinFillN(rdone_,nrows1,PRESENT_IN_REDUCED) ;
1390  CoinZeroN(rdone_+nrows1,nrows0_in-nrows1) ;
1391# endif
1392}
1393/* This is main part of Presolve */
1394ClpSimplex * 
1395ClpPresolve::gutsOfPresolvedModel(ClpSimplex * originalModel,
1396                                  double feasibilityTolerance,
1397                                  bool keepIntegers,
1398                                  int numberPasses,
1399                                  bool dropNames,
1400                                  bool doRowObjective)
1401{
1402  ncols_ = originalModel->getNumCols();
1403  nrows_ = originalModel->getNumRows();
1404  nelems_ = originalModel->getNumElements();
1405  numberPasses_ = numberPasses;
1406
1407  double maxmin = originalModel->getObjSense();
1408  originalModel_ = originalModel;
1409  delete [] originalColumn_;
1410  originalColumn_ = new int[ncols_];
1411  delete [] originalRow_;
1412  originalRow_ = new int[nrows_];
1413  // and fill in case returns early
1414  int i;
1415  for (i=0;i<ncols_;i++)
1416    originalColumn_[i]=i;
1417  for (i=0;i<nrows_;i++)
1418    originalRow_[i]=i;
1419  delete [] rowObjective_;
1420  if (doRowObjective) {
1421    rowObjective_ = new double [nrows_];
1422    memset(rowObjective_,0,nrows_*sizeof(double));
1423  } else {
1424    rowObjective_=NULL;
1425  }
1426
1427  // result is 0 - okay, 1 infeasible, -1 go round again, 2 - original model
1428  int result = -1;
1429 
1430  // User may have deleted - its their responsibility
1431  presolvedModel_=NULL;
1432  // Messages
1433  CoinMessages messages = originalModel->coinMessages();
1434  // Only go round 100 times even if integer preprocessing
1435  int totalPasses=100;
1436  while (result==-1) {
1437
1438#ifndef CLP_NO_STD
1439    // make new copy
1440    if (saveFile_=="") {
1441#endif
1442      delete presolvedModel_;
1443#ifndef CLP_NO_STD
1444      // So won't get names
1445      int lengthNames = originalModel->lengthNames();
1446      originalModel->setLengthNames(0);
1447#endif
1448      presolvedModel_ = new ClpSimplex(*originalModel);
1449#ifndef CLP_NO_STD
1450      originalModel->setLengthNames(lengthNames);
1451    } else {
1452      presolvedModel_=originalModel;
1453    }
1454    presolvedModel_->dropNames();
1455#endif
1456
1457    // drop integer information if wanted
1458    if (!keepIntegers)
1459      presolvedModel_->deleteIntegerInformation();
1460    totalPasses--;
1461
1462    double ratio=2.0;
1463    if (substitution_>3)
1464      ratio=substitution_;
1465    else if (substitution_==2)
1466      ratio=1.5;
1467    CoinPresolveMatrix prob(ncols_,
1468                        maxmin,
1469                        presolvedModel_,
1470                        nrows_, nelems_,true,nonLinearValue_,ratio);
1471    prob.setMaximumSubstitutionLevel(substitution_);
1472    if (doRowObjective) 
1473      memset(rowObjective_,0,nrows_*sizeof(double));
1474    // See if we want statistics
1475    if ((presolveActions_&0x80000000)!=0)
1476      prob.statistics();
1477    // make sure row solution correct
1478    {
1479      double *colels    = prob.colels_;
1480      int *hrow         = prob.hrow_;
1481      CoinBigIndex *mcstrt              = prob.mcstrt_;
1482      int *hincol               = prob.hincol_;
1483      int ncols         = prob.ncols_;
1484     
1485     
1486      double * csol = prob.sol_;
1487      double * acts = prob.acts_;
1488      int nrows = prob.nrows_;
1489
1490      int colx;
1491
1492      memset(acts,0,nrows*sizeof(double));
1493     
1494      for (colx = 0; colx < ncols; ++colx) {
1495        double solutionValue = csol[colx];
1496        for (int i=mcstrt[colx]; i<mcstrt[colx]+hincol[colx]; ++i) {
1497          int row = hrow[i];
1498          double coeff = colels[i];
1499          acts[row] += solutionValue*coeff;
1500        }
1501      }
1502    }
1503
1504    // move across feasibility tolerance
1505    prob.feasibilityTolerance_ = feasibilityTolerance;
1506
1507    // Do presolve
1508    paction_ = presolve(&prob);
1509
1510    result =0; 
1511
1512    if (prob.status_==0&&paction_) {
1513      // Looks feasible but double check to see if anything slipped through
1514      int n             = prob.ncols_;
1515      double * lo = prob.clo_;
1516      double * up = prob.cup_;
1517      int i;
1518     
1519      for (i=0;i<n;i++) {
1520        if (up[i]<lo[i]) {
1521          if (up[i]<lo[i]-1.0e-8) {
1522            // infeasible
1523            prob.status_=1;
1524          } else {
1525            up[i]=lo[i];
1526          }
1527        }
1528      }
1529     
1530      n = prob.nrows_;
1531      lo = prob.rlo_;
1532      up = prob.rup_;
1533
1534      for (i=0;i<n;i++) {
1535        if (up[i]<lo[i]) {
1536          if (up[i]<lo[i]-1.0e-8) {
1537            // infeasible
1538            prob.status_=1;
1539          } else {
1540            up[i]=lo[i];
1541          }
1542        }
1543      }
1544    }
1545    if (prob.status_==0&&paction_) {
1546      // feasible
1547   
1548      prob.update_model(presolvedModel_, nrows_, ncols_, nelems_);
1549      // copy status and solution
1550      CoinMemcpyN(           prob.sol_,prob.ncols_,presolvedModel_->primalColumnSolution());
1551      CoinMemcpyN(           prob.acts_,prob.nrows_,presolvedModel_->primalRowSolution());
1552      CoinMemcpyN(           prob.colstat_,prob.ncols_,presolvedModel_->statusArray());
1553      CoinMemcpyN(           prob.rowstat_,prob.nrows_,presolvedModel_->statusArray()+prob.ncols_);
1554      delete [] prob.sol_;
1555      delete [] prob.acts_;
1556      delete [] prob.colstat_;
1557      prob.sol_=NULL;
1558      prob.acts_=NULL;
1559      prob.colstat_=NULL;
1560     
1561      int ncolsNow = presolvedModel_->getNumCols();
1562      CoinMemcpyN(prob.originalColumn_,ncolsNow,originalColumn_);
1563#ifndef SLIM_CLP
1564#ifndef NO_RTTI
1565      ClpQuadraticObjective * quadraticObj = (dynamic_cast< ClpQuadraticObjective*>(originalModel->objectiveAsObject()));
1566#else
1567      ClpQuadraticObjective * quadraticObj = NULL;
1568      if (originalModel->objectiveAsObject()->type()==2)
1569        quadraticObj = (static_cast< ClpQuadraticObjective*>(originalModel->objectiveAsObject()));
1570#endif
1571      if (quadraticObj) {
1572        // set up for subset
1573        char * mark = new char [ncols_];
1574        memset(mark,0,ncols_);
1575        CoinPackedMatrix * quadratic = quadraticObj->quadraticObjective();
1576        //const int * columnQuadratic = quadratic->getIndices();
1577        //const CoinBigIndex * columnQuadraticStart = quadratic->getVectorStarts();
1578        const int * columnQuadraticLength = quadratic->getVectorLengths();
1579        //double * quadraticElement = quadratic->getMutableElements();
1580        int numberColumns = quadratic->getNumCols();
1581        ClpQuadraticObjective * newObj = new ClpQuadraticObjective(*quadraticObj,
1582                                                                   ncolsNow,
1583                                                                   originalColumn_);
1584        // and modify linear and check
1585        double * linear = newObj->linearObjective();
1586 CoinMemcpyN(presolvedModel_->objective(),ncolsNow,linear);
1587        int iColumn;
1588        for ( iColumn=0;iColumn<numberColumns;iColumn++) {
1589          if (columnQuadraticLength[iColumn])
1590            mark[iColumn]=1;
1591        }
1592        // and new
1593        quadratic = newObj->quadraticObjective();
1594        columnQuadraticLength = quadratic->getVectorLengths();
1595        int numberColumns2 = quadratic->getNumCols();
1596        for ( iColumn=0;iColumn<numberColumns2;iColumn++) {
1597          if (columnQuadraticLength[iColumn])
1598            mark[originalColumn_[iColumn]]=0;
1599        }
1600        presolvedModel_->setObjective(newObj);
1601        delete newObj;
1602        // final check
1603        for ( iColumn=0;iColumn<numberColumns;iColumn++) 
1604          if (mark[iColumn])
1605            printf("Quadratic column %d modified - may be okay\n",iColumn);
1606        delete [] mark;
1607      }
1608#endif
1609      delete [] prob.originalColumn_;
1610      prob.originalColumn_=NULL;
1611      int nrowsNow = presolvedModel_->getNumRows();
1612      CoinMemcpyN(prob.originalRow_,nrowsNow,originalRow_);
1613      delete [] prob.originalRow_;
1614      prob.originalRow_=NULL;
1615#ifndef CLP_NO_STD
1616      if (!dropNames&&originalModel->lengthNames()) {
1617        // Redo names
1618        int iRow;
1619        std::vector<std::string> rowNames;
1620        rowNames.reserve(nrowsNow);
1621        for (iRow=0;iRow<nrowsNow;iRow++) {
1622          int kRow = originalRow_[iRow];
1623          rowNames.push_back(originalModel->rowName(kRow));
1624        }
1625     
1626        int iColumn;
1627        std::vector<std::string> columnNames;
1628        columnNames.reserve(ncolsNow);
1629        for (iColumn=0;iColumn<ncolsNow;iColumn++) {
1630          int kColumn = originalColumn_[iColumn];
1631          columnNames.push_back(originalModel->columnName(kColumn));
1632        }
1633        presolvedModel_->copyNames(rowNames,columnNames);
1634      } else {
1635        presolvedModel_->setLengthNames(0);
1636      }
1637#endif
1638      if (rowObjective_) {
1639        int iRow;
1640        int k=-1;
1641        int nObj=0;
1642        for (iRow=0;iRow<nrowsNow;iRow++) {
1643          int kRow = originalRow_[iRow];
1644          assert (kRow>k);
1645          k=kRow;
1646          rowObjective_[iRow]=rowObjective_[kRow];
1647          if (rowObjective_[iRow])
1648            nObj++;
1649        }
1650        if (nObj) {
1651          printf("%d costed slacks\n",nObj);
1652          presolvedModel_->setRowObjective(rowObjective_);
1653        }
1654      }
1655      // now clean up integer variables.  This can modify original
1656      int i;
1657      const char * information = presolvedModel_->integerInformation();
1658      if (information) {
1659        int numberChanges=0;
1660        double * lower0 = originalModel_->columnLower();
1661        double * upper0 = originalModel_->columnUpper();
1662        double * lower = presolvedModel_->columnLower();
1663        double * upper = presolvedModel_->columnUpper();
1664        for (i=0;i<ncolsNow;i++) {
1665          if (!information[i])
1666            continue;
1667          int iOriginal = originalColumn_[i];
1668          double lowerValue0 = lower0[iOriginal];
1669          double upperValue0 = upper0[iOriginal];
1670          double lowerValue = ceil(lower[i]-1.0e-5);
1671          double upperValue = floor(upper[i]+1.0e-5);
1672          lower[i]=lowerValue;
1673          upper[i]=upperValue;
1674          if (lowerValue>upperValue) {
1675            numberChanges++;
1676            presolvedModel_->messageHandler()->message(COIN_PRESOLVE_COLINFEAS,
1677                                                       messages)
1678                                                         <<iOriginal
1679                                                         <<lowerValue
1680                                                         <<upperValue
1681                                                         <<CoinMessageEol;
1682            result=1;
1683          } else {
1684            if (lowerValue>lowerValue0+1.0e-8) {
1685              lower0[iOriginal] = lowerValue;
1686              numberChanges++;
1687            }
1688            if (upperValue<upperValue0-1.0e-8) {
1689              upper0[iOriginal] = upperValue;
1690              numberChanges++;
1691            }
1692          }       
1693        }
1694        if (numberChanges) {
1695          presolvedModel_->messageHandler()->message(COIN_PRESOLVE_INTEGERMODS,
1696                                                     messages)
1697                                                       <<numberChanges
1698                                                       <<CoinMessageEol;
1699          if (!result&&totalPasses>0) {
1700            result = -1; // round again
1701            const CoinPresolveAction *paction = paction_;
1702            while (paction) {
1703              const CoinPresolveAction *next = paction->next;
1704              delete paction;
1705              paction = next;
1706            }
1707            paction_=NULL;
1708          }
1709        }
1710      }
1711    } else if (prob.status_) {
1712      // infeasible or unbounded
1713      result=1;
1714      originalModel->setProblemStatus(prob.status_);
1715    } else {
1716      // no changes - model needs restoring after Lou's changes
1717#ifndef CLP_NO_STD
1718      if (saveFile_=="") {
1719#endif
1720        delete presolvedModel_;
1721        presolvedModel_ = new ClpSimplex(*originalModel);
1722        // but we need to remove gaps
1723        ClpPackedMatrix* clpMatrix =
1724          dynamic_cast< ClpPackedMatrix*>(presolvedModel_->clpMatrix());
1725        if (clpMatrix) {
1726          clpMatrix->getPackedMatrix()->removeGaps();
1727        }
1728#ifndef CLP_NO_STD
1729      } else {
1730        presolvedModel_=originalModel;
1731      }
1732      presolvedModel_->dropNames();
1733#endif
1734     
1735      // drop integer information if wanted
1736      if (!keepIntegers)
1737        presolvedModel_->deleteIntegerInformation();
1738      result=2;
1739    }
1740  }
1741  if (result==0||result==2) {
1742    int nrowsAfter = presolvedModel_->getNumRows();
1743    int ncolsAfter = presolvedModel_->getNumCols();
1744    CoinBigIndex nelsAfter = presolvedModel_->getNumElements();
1745    presolvedModel_->messageHandler()->message(COIN_PRESOLVE_STATS,
1746                                               messages)
1747                                                 <<nrowsAfter<< -(nrows_ - nrowsAfter)
1748                                                 << ncolsAfter<< -(ncols_ - ncolsAfter)
1749                                                 <<nelsAfter<< -(nelems_ - nelsAfter)
1750                                                 <<CoinMessageEol;
1751  } else {
1752    destroyPresolve();
1753    if (presolvedModel_!=originalModel_)
1754      delete presolvedModel_;
1755    presolvedModel_=NULL;
1756  }
1757  return presolvedModel_;
1758}
1759
1760
Note: See TracBrowser for help on using the repository browser.