source: trunk/ClpPresolve.cpp @ 701

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

messages and preprocessing

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