source: trunk/ClpPresolve.cpp @ 461

Last change on this file since 461 was 461, checked in by forrest, 15 years ago

Trying to make faster

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 42.0 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#include "ClpQuadraticObjective.hpp"
17
18#include "ClpPresolve.hpp"
19#include "CoinPresolveMatrix.hpp"
20
21#include "CoinPresolveEmpty.hpp"
22#include "CoinPresolveFixed.hpp"
23#include "CoinPresolvePsdebug.hpp"
24#include "CoinPresolveSingleton.hpp"
25#include "CoinPresolveDoubleton.hpp"
26#include "CoinPresolveTripleton.hpp"
27#include "CoinPresolveZeros.hpp"
28#include "CoinPresolveSubst.hpp"
29#include "CoinPresolveForcing.hpp"
30#include "CoinPresolveDual.hpp"
31#include "CoinPresolveTighten.hpp"
32#include "CoinPresolveUseless.hpp"
33#include "CoinPresolveDupcol.hpp"
34#include "CoinPresolveImpliedFree.hpp"
35#include "CoinPresolveIsolated.hpp"
36#include "CoinMessage.hpp"
37
38
39
40ClpPresolve::ClpPresolve() :
41  originalModel_(NULL),
42  presolvedModel_(NULL),
43  nonLinearValue_(0.0),
44  originalColumn_(NULL),
45  originalRow_(NULL),
46  paction_(0),
47  ncols_(0),
48  nrows_(0),
49  nelems_(0),
50  numberPasses_(5),
51  saveFile_(""),
52  presolveActions_(0)
53{
54}
55
56ClpPresolve::~ClpPresolve()
57{
58  destroyPresolve();
59}
60// Gets rid of presolve actions (e.g.when infeasible)
61void 
62ClpPresolve::destroyPresolve()
63{
64 const CoinPresolveAction *paction = paction_;
65  while (paction) {
66    const CoinPresolveAction *next = paction->next;
67    delete paction;
68    paction = next;
69  }
70  delete [] originalColumn_;
71  delete [] originalRow_;
72  paction_=NULL;
73  originalColumn_=NULL;
74  originalRow_=NULL;
75}
76
77/* This version of presolve returns a pointer to a new presolved
78   model.  NULL if infeasible
79*/
80ClpSimplex * 
81ClpPresolve::presolvedModel(ClpSimplex & si,
82                         double feasibilityTolerance,
83                         bool keepIntegers,
84                         int numberPasses,
85                         bool dropNames)
86{
87  // Check matrix
88  if (!si.clpMatrix()->allElementsInRange(&si,si.getSmallElementValue(),
89                                          1.0e20))
90    return NULL;
91  else
92    return gutsOfPresolvedModel(&si,feasibilityTolerance,keepIntegers,numberPasses,dropNames);
93}
94/* This version of presolve updates
95   model and saves original data to file.  Returns non-zero if infeasible
96*/
97int
98ClpPresolve::presolvedModelToFile(ClpSimplex &si,std::string fileName,
99                            double feasibilityTolerance,
100                            bool keepIntegers,
101                            int numberPasses)
102{
103  // Check matrix
104  if (!si.clpMatrix()->allElementsInRange(&si,si.getSmallElementValue(),
105                                          1.0e20))
106    return 2;
107  saveFile_=fileName;
108  si.saveModel(saveFile_.c_str());
109  ClpSimplex * model = gutsOfPresolvedModel(&si,feasibilityTolerance,keepIntegers,numberPasses,true);
110  if (model==&si) {
111    return 0;
112  } else {
113    si.restoreModel(saveFile_.c_str());
114    return 1;
115  }
116}
117
118// Return pointer to presolved model
119ClpSimplex * 
120ClpPresolve::model() const
121{
122  return presolvedModel_;
123}
124// Return pointer to original model
125ClpSimplex * 
126ClpPresolve::originalModel() const
127{
128  return originalModel_;
129}
130void 
131ClpPresolve::postsolve(bool updateStatus)
132{
133  // Return at once if no presolved model
134  if (!presolvedModel_)
135    return;
136  // Messages
137  CoinMessages messages = CoinMessage(presolvedModel_->messages().language());
138  if (!presolvedModel_->isProvenOptimal()) {
139    presolvedModel_->messageHandler()->message(COIN_PRESOLVE_NONOPTIMAL,
140                                             messages)
141                                               <<CoinMessageEol;
142  }
143
144  // this is the size of the original problem
145  const int ncols0  = ncols_;
146  const int nrows0  = nrows_;
147  const CoinBigIndex nelems0 = nelems_;
148 
149  // this is the reduced problem
150  int ncols = presolvedModel_->getNumCols();
151  int nrows = presolvedModel_->getNumRows();
152
153  double * acts=NULL;
154  double * sol =NULL;
155  unsigned char * rowstat=NULL;
156  unsigned char * colstat = NULL;
157  if (saveFile_=="") {
158    // reality check
159    assert(ncols0==originalModel_->getNumCols());
160    assert(nrows0==originalModel_->getNumRows());
161    acts = originalModel_->primalRowSolution();
162    sol  = originalModel_->primalColumnSolution();
163    if (updateStatus) {
164      unsigned char *status = originalModel_->statusArray();
165      rowstat = status + ncols0;
166      colstat = status;
167      memcpy(colstat, presolvedModel_->statusArray(), ncols);
168      memcpy(rowstat, presolvedModel_->statusArray()+ncols, nrows);
169    }
170  } else {
171    // from file
172    acts = new double[nrows0];
173    sol  = new double[ncols0];
174    CoinZeroN(acts,nrows0);
175    CoinZeroN(sol,ncols0);
176    if (updateStatus) {
177      unsigned char *status = new unsigned char [nrows0+ncols0];
178      rowstat = status + ncols0;
179      colstat = status;
180      memcpy(colstat, presolvedModel_->statusArray(), ncols);
181      memcpy(rowstat, presolvedModel_->statusArray()+ncols, nrows);
182    }
183  }
184
185  // CoinPostsolveMatrix object assumes ownership of sol, acts, colstat;
186  // will be deleted by ~CoinPostsolveMatrix. delete[] operations below
187  // cause duplicate free. In case where saveFile == "", as best I can see
188  // arrays are owned by originalModel_. fix is to
189  // clear fields in prob to avoid delete[] in ~CoinPostsolveMatrix.
190  CoinPostsolveMatrix prob(presolvedModel_,
191                       ncols0,
192                       nrows0,
193                       nelems0,
194                       presolvedModel_->getObjSense(),
195                       // end prepost
196                       
197                       sol, acts,
198                       colstat, rowstat);
199   
200  postsolve(prob);
201
202  if (saveFile_!="") {
203    // From file
204    assert (originalModel_==presolvedModel_);
205    originalModel_->restoreModel(saveFile_.c_str());
206    memcpy(originalModel_->primalRowSolution(),acts,nrows0*sizeof(double));
207    // delete [] acts;
208    memcpy(originalModel_->primalColumnSolution(),sol,ncols0*sizeof(double));
209    // delete [] sol;
210    if (updateStatus) {
211      memcpy(originalModel_->statusArray(),colstat,nrows0+ncols0);
212      // delete [] colstat;
213    }
214  } else {
215    prob.sol_ = 0 ;
216    prob.acts_ = 0 ;
217    prob.colstat_ = 0 ;
218  }
219  // put back duals
220  memcpy(originalModel_->dualRowSolution(),prob.rowduals_,
221         nrows_*sizeof(double));
222  double maxmin = originalModel_->getObjSense();
223  if (maxmin<0.0) {
224    // swap signs
225    int i;
226    double * pi = originalModel_->dualRowSolution();
227    for (i=0;i<nrows_;i++)
228      pi[i] = -pi[i];
229  }
230  // Now check solution
231  double offset;
232  memcpy(originalModel_->dualColumnSolution(),
233         originalModel_->objectiveAsObject()->gradient(originalModel_,
234                                                       originalModel_->primalColumnSolution(),
235                                                       offset,true),ncols_*sizeof(double));
236  originalModel_->transposeTimes(-1.0,
237                                 originalModel_->dualRowSolution(),
238                                 originalModel_->dualColumnSolution());
239  memset(originalModel_->primalRowSolution(),0,nrows_*sizeof(double));
240  originalModel_->times(1.0,originalModel_->primalColumnSolution(),
241                        originalModel_->primalRowSolution());
242  originalModel_->checkSolution();
243  // Messages
244  presolvedModel_->messageHandler()->message(COIN_PRESOLVE_POSTSOLVE,
245                                            messages)
246                                              <<originalModel_->objectiveValue()
247                                              <<originalModel_->sumDualInfeasibilities()
248                                              <<originalModel_->numberDualInfeasibilities()
249                                              <<originalModel_->sumPrimalInfeasibilities()
250                                              <<originalModel_->numberPrimalInfeasibilities()
251                                               <<CoinMessageEol;
252 
253  //originalModel_->objectiveValue_=objectiveValue_;
254  originalModel_->setNumberIterations(presolvedModel_->numberIterations());
255  if (!presolvedModel_->status()) {
256    if (!originalModel_->numberDualInfeasibilities()&&
257        !originalModel_->numberPrimalInfeasibilities()) {
258      originalModel_->setProblemStatus( 0);
259    } else {
260      originalModel_->setProblemStatus( -1);
261      presolvedModel_->messageHandler()->message(COIN_PRESOLVE_NEEDS_CLEANING,
262                                            messages)
263                                              <<CoinMessageEol;
264    }
265  } else {
266    originalModel_->setProblemStatus( presolvedModel_->status());
267  }
268  if (saveFile_!="") 
269    presolvedModel_=NULL;
270}
271
272// return pointer to original columns
273const int * 
274ClpPresolve::originalColumns() const
275{
276  return originalColumn_;
277}
278// return pointer to original rows
279const int * 
280ClpPresolve::originalRows() const
281{
282  return originalRow_;
283}
284// Set pointer to original model
285void 
286ClpPresolve::setOriginalModel(ClpSimplex * model)
287{
288  originalModel_=model;
289}
290#if 0
291// A lazy way to restrict which transformations are applied
292// during debugging.
293static int ATOI(const char *name)
294{
295 return true;
296#if     PRESOLVE_DEBUG || PRESOLVE_SUMMARY
297  if (getenv(name)) {
298    int val = atoi(getenv(name));
299    printf("%s = %d\n", name, val);
300    return (val);
301  } else {
302    if (strcmp(name,"off"))
303      return (true);
304    else
305      return (false);
306  }
307#else
308  return (true);
309#endif
310}
311#endif
312//#define PRESOLVE_DEBUG 1
313#if PRESOLVE_DEBUG
314void check_sol(CoinPresolveMatrix *prob,double tol)
315{
316  double *colels        = prob->colels_;
317  int *hrow             = prob->hrow_;
318  int *mcstrt           = prob->mcstrt_;
319  int *hincol           = prob->hincol_;
320  int *hinrow           = prob->hinrow_;
321  int ncols             = prob->ncols_;
322
323
324  double * csol = prob->sol_;
325  double * acts = prob->acts_;
326  double * clo = prob->clo_;
327  double * cup = prob->cup_;
328  int nrows = prob->nrows_;
329  double * rlo = prob->rlo_;
330  double * rup = prob->rup_;
331
332  int colx;
333
334  double * rsol = new double[nrows];
335  memset(rsol,0,nrows*sizeof(double));
336
337  for (colx = 0; colx < ncols; ++colx) {
338    if (1) {
339      CoinBigIndex k = mcstrt[colx];
340      int nx = hincol[colx];
341      double solutionValue = csol[colx];
342      for (int i=0; i<nx; ++i) {
343        int row = hrow[k];
344        double coeff = colels[k];
345        k++;
346        rsol[row] += solutionValue*coeff;
347      }
348      if (csol[colx]<clo[colx]-tol) {
349        printf("low CSOL:  %d  - %g %g %g\n",
350                   colx, clo[colx], csol[colx], cup[colx]);
351      } else if (csol[colx]>cup[colx]+tol) {
352        printf("high CSOL:  %d  - %g %g %g\n",
353                   colx, clo[colx], csol[colx], cup[colx]);
354      } 
355    }
356  }
357  int rowx;
358  for (rowx = 0; rowx < nrows; ++rowx) {
359    if (hinrow[rowx]) {
360      if (fabs(rsol[rowx]-acts[rowx])>tol)
361        printf("inacc RSOL:  %d - %g %g (acts_ %g) %g\n",
362                   rowx,  rlo[rowx], rsol[rowx], acts[rowx], rup[rowx]);
363      if (rsol[rowx]<rlo[rowx]-tol) {
364        printf("low RSOL:  %d - %g %g %g\n",
365                   rowx,  rlo[rowx], rsol[rowx], rup[rowx]);
366      } else if (rsol[rowx]>rup[rowx]+tol ) {
367        printf("high RSOL:  %d - %g %g %g\n",
368                   rowx,  rlo[rowx], rsol[rowx], rup[rowx]);
369      } 
370    }
371  }
372  delete [] rsol;
373}
374#endif
375// This is the presolve loop.
376// It is a separate virtual function so that it can be easily
377// customized by subclassing CoinPresolve.
378const CoinPresolveAction *ClpPresolve::presolve(CoinPresolveMatrix *prob)
379{
380  // Messages
381  CoinMessages messages = CoinMessage(prob->messages().language());
382  paction_ = 0;
383
384  prob->status_=0; // say feasible
385
386  paction_ = make_fixed(prob, paction_);
387  // if integers then switch off dual stuff
388  // later just do individually
389  bool doDualStuff = (presolvedModel_->integerInformation()==NULL);
390  if (prob->anyProhibited()) 
391    doDualStuff=false;
392  if (!doDual())
393    doDualStuff=false;
394#if     PRESOLVE_CONSISTENCY
395//  presolve_links_ok(prob->rlink_, prob->mrstrt_, prob->hinrow_, prob->nrows_);
396    presolve_links_ok(prob,false,true) ;
397#endif
398
399  if (!prob->status_) {
400    const bool slackd = doSingleton();
401    const bool doubleton = doDoubleton();
402    const bool tripleton = doTripleton();
403    const bool forcing = doForcing();
404    const bool ifree = doImpliedFree();
405    const bool zerocost = doTighten();
406    const bool dupcol = doDupcol();
407    const bool duprow = doDuprow();
408    const bool dual = doDualStuff;
409   
410    // some things are expensive so just do once (normally)
411
412    int i;
413    // say look at all
414    if (!prob->anyProhibited()) {
415      for (i=0;i<nrows_;i++) 
416        prob->rowsToDo_[i]=i;
417      prob->numberRowsToDo_=nrows_;
418      for (i=0;i<ncols_;i++) 
419        prob->colsToDo_[i]=i;
420      prob->numberColsToDo_=ncols_;
421    } else {
422      // some stuff must be left alone
423      prob->numberRowsToDo_=0;
424      for (i=0;i<nrows_;i++) 
425        if (!prob->rowProhibited(i))
426            prob->rowsToDo_[prob->numberRowsToDo_++]=i;
427      prob->numberColsToDo_=0;
428      for (i=0;i<ncols_;i++) 
429        if (!prob->colProhibited(i))
430            prob->colsToDo_[prob->numberColsToDo_++]=i;
431    }
432
433
434    int iLoop;
435#if     PRESOLVE_DEBUG
436    check_sol(prob,1.0e0);
437#endif
438    if (dupcol) {
439      //paction_ = dupcol_action::presolve(prob, paction_);
440    }
441
442    // Check number rows dropped
443    int lastDropped=0;
444    prob->pass_=0;
445    for (iLoop=0;iLoop<numberPasses_;iLoop++) {
446    // See if we want statistics
447    if ((presolveActions_&0x80000000)!=0)
448      printf("Starting major pass %d after %g seconds\n",iLoop+1,CoinCpuTime()-prob->startTime_);
449#ifdef PRESOLVE_SUMMARY
450      printf("Starting major pass %d\n",iLoop+1);
451#endif
452      const CoinPresolveAction * const paction0 = paction_;
453      // look for substitutions with no fill
454      int fill_level=2;
455      //fill_level=10;
456      //printf("** fill_level == 10 !\n");
457      int whichPass=0;
458      while (1) {
459        whichPass++;
460        prob->pass_++;
461        const CoinPresolveAction * const paction1 = paction_;
462
463        if (slackd) {
464          bool notFinished = true;
465          while (notFinished) 
466            paction_ = slack_doubleton_action::presolve(prob, paction_,
467                                                        notFinished);
468          if (prob->status_)
469            break;
470        }
471        if (dual&&whichPass==1) {
472          // this can also make E rows so do one bit here
473          paction_ = remove_dual_action::presolve(prob, paction_);
474          if (prob->status_)
475            break;
476        }
477
478        if (doubleton) {
479          paction_ = doubleton_action::presolve(prob, paction_);
480          if (prob->status_)
481            break;
482        }
483
484        if (tripleton) {
485          paction_ = tripleton_action::presolve(prob, paction_);
486          if (prob->status_)
487            break;
488        }
489
490        if (zerocost) {
491          paction_ = do_tighten_action::presolve(prob, paction_);
492          if (prob->status_)
493            break;
494        }
495
496        if (forcing) {
497          paction_ = forcing_constraint_action::presolve(prob, paction_);
498          if (prob->status_)
499            break;
500        }
501
502        if (ifree&&(whichPass%5)==1) {
503        paction_ = implied_free_action::presolve(prob, paction_,fill_level);
504        if (prob->status_)
505          break;
506        }
507
508#if     PRESOLVE_DEBUG
509        check_sol(prob,1.0e0);
510#endif
511
512#if     PRESOLVE_CONSISTENCY
513//      presolve_links_ok(prob->rlink_, prob->mrstrt_, prob->hinrow_,
514//                        prob->nrows_);
515        presolve_links_ok(prob,false,true) ;
516#endif
517
518//#if   PRESOLVE_DEBUG
519//      presolve_no_zeros(prob->mcstrt_, prob->colels_, prob->hincol_,
520//                        prob->ncols_);
521//#endif
522//#if   PRESOLVE_CONSISTENCY
523//      prob->consistent();
524//#endif
525#if     PRESOLVE_CONSISTENCY
526        presolve_no_zeros(prob,true,false) ;
527        presolve_consistent(prob,true) ;
528#endif
529
530         
531        // set up for next pass
532        // later do faster if many changes i.e. memset and memcpy
533        prob->numberRowsToDo_ = prob->numberNextRowsToDo_;
534        int kcheck;
535        bool found=false;
536        kcheck=-1;
537        for (i=0;i<prob->numberNextRowsToDo_;i++) {
538          int index = prob->nextRowsToDo_[i];
539          prob->unsetRowChanged(index);
540          prob->rowsToDo_[i] = index;
541          if (index==kcheck) {
542            printf("row %d on list after pass %d\n",kcheck,
543                   whichPass);
544            found=true;
545          }
546        }
547        if (!found&&kcheck>=0)
548          prob->rowsToDo_[prob->numberRowsToDo_++]=kcheck;
549        prob->numberNextRowsToDo_=0;
550        prob->numberColsToDo_ = prob->numberNextColsToDo_;
551        kcheck=-1;
552        found=false;
553        for (i=0;i<prob->numberNextColsToDo_;i++) {
554          int index = prob->nextColsToDo_[i];
555          prob->unsetColChanged(index);
556          prob->colsToDo_[i] = index;
557          if (index==kcheck) {
558            printf("col %d on list after pass %d\n",kcheck,
559                   whichPass);
560            found=true;
561          }
562        }
563        if (!found&&kcheck>=0)
564          prob->colsToDo_[prob->numberColsToDo_++]=kcheck;
565        prob->numberNextColsToDo_=0;
566        if (paction_ == paction1&&fill_level>0)
567          break;
568      }
569      // say look at all
570      int i;
571      if (!prob->anyProhibited()) {
572        for (i=0;i<nrows_;i++) 
573          prob->rowsToDo_[i]=i;
574        prob->numberRowsToDo_=nrows_;
575        for (i=0;i<ncols_;i++) 
576          prob->colsToDo_[i]=i;
577        prob->numberColsToDo_=ncols_;
578      } else {
579        // some stuff must be left alone
580        prob->numberRowsToDo_=0;
581        for (i=0;i<nrows_;i++) 
582          if (!prob->rowProhibited(i))
583            prob->rowsToDo_[prob->numberRowsToDo_++]=i;
584        prob->numberColsToDo_=0;
585        for (i=0;i<ncols_;i++) 
586          if (!prob->colProhibited(i))
587            prob->colsToDo_[prob->numberColsToDo_++]=i;
588      }
589      // now expensive things
590      // this caused world.mps to run into numerical difficulties
591#ifdef PRESOLVE_SUMMARY
592      printf("Starting expensive\n");
593#endif
594
595      if (dual) {
596        int itry;
597        for (itry=0;itry<5;itry++) {
598          paction_ = remove_dual_action::presolve(prob, paction_);
599          if (prob->status_)
600            break;
601          const CoinPresolveAction * const paction2 = paction_;
602          if (ifree) {
603            //int fill_level=0; // switches off substitution
604            paction_ = implied_free_action::presolve(prob, paction_,fill_level);
605            if (prob->status_)
606              break;
607          }
608          if (paction_ == paction2)
609            break;
610        }
611      } else if (ifree) {
612        // just free
613        int fill_level=0; // switches off substitution
614        paction_ = implied_free_action::presolve(prob, paction_,fill_level);
615        if (prob->status_)
616          break;
617      }
618#if     PRESOLVE_DEBUG
619      check_sol(prob,1.0e0);
620#endif
621      if (dupcol) {
622        paction_ = dupcol_action::presolve(prob, paction_);
623        if (prob->status_)
624          break;
625      }
626#if     PRESOLVE_DEBUG
627        check_sol(prob,1.0e0);
628#endif
629     
630      if (duprow) {
631        paction_ = duprow_action::presolve(prob, paction_);
632        if (prob->status_)
633          break;
634      }
635#if     PRESOLVE_DEBUG
636      check_sol(prob,1.0e0);
637#endif
638      {
639        int * hinrow = prob->hinrow_;
640        int numberDropped=0;
641        for (i=0;i<nrows_;i++) 
642          if (!hinrow[i])
643            numberDropped++;
644       
645        prob->messageHandler()->message(COIN_PRESOLVE_PASS,
646                                        messages)
647                                          <<numberDropped<<iLoop+1
648                                          <<CoinMessageEol;
649        //printf("%d rows dropped after pass %d\n",numberDropped,
650        //     iLoop+1);
651        if (numberDropped==lastDropped)
652          break;
653        else
654          lastDropped = numberDropped;
655      }
656      if (paction_ == paction0)
657        break;
658         
659    }
660  }
661  if (!prob->status_) {
662    paction_ = drop_zero_coefficients(prob, paction_);
663#if     PRESOLVE_DEBUG
664        check_sol(prob,1.0e0);
665#endif
666
667    paction_ = drop_empty_cols_action::presolve(prob, paction_);
668    paction_ = drop_empty_rows_action::presolve(prob, paction_);
669#if     PRESOLVE_DEBUG
670        check_sol(prob,1.0e0);
671#endif
672  }
673 
674  if (prob->status_) {
675    if (prob->status_==1)
676          prob->messageHandler()->message(COIN_PRESOLVE_INFEAS,
677                                             messages)
678                                               <<prob->feasibilityTolerance_
679                                               <<CoinMessageEol;
680    else if (prob->status_==2)
681          prob->messageHandler()->message(COIN_PRESOLVE_UNBOUND,
682                                             messages)
683                                               <<CoinMessageEol;
684    else
685          prob->messageHandler()->message(COIN_PRESOLVE_INFEASUNBOUND,
686                                             messages)
687                                               <<CoinMessageEol;
688    // get rid of data
689    destroyPresolve();
690  }
691  return (paction_);
692}
693
694void check_djs(CoinPostsolveMatrix *prob);
695
696
697// We could have implemented this by having each postsolve routine
698// directly call the next one, but this may make it easier to add debugging checks.
699void ClpPresolve::postsolve(CoinPostsolveMatrix &prob)
700{
701  {
702    // Check activities
703    double *colels      = prob.colels_;
704    int *hrow           = prob.hrow_;
705    CoinBigIndex *mcstrt                = prob.mcstrt_;
706    int *hincol         = prob.hincol_;
707    int *link           = prob.link_;
708    int ncols           = prob.ncols_;
709
710    char *cdone = prob.cdone_;
711
712    double * csol = prob.sol_;
713    int nrows = prob.nrows_;
714
715    int colx;
716
717    double * rsol = prob.acts_;
718    memset(rsol,0,nrows*sizeof(double));
719
720    for (colx = 0; colx < ncols; ++colx) {
721      if (cdone[colx]) {
722        CoinBigIndex k = mcstrt[colx];
723        int nx = hincol[colx];
724        double solutionValue = csol[colx];
725        for (int i=0; i<nx; ++i) {
726          int row = hrow[k];
727          double coeff = colels[k];
728          k = link[k];
729          rsol[row] += solutionValue*coeff;
730        }
731      }
732    }
733  }
734  const CoinPresolveAction *paction = paction_;
735
736#if     PRESOLVE_DEBUG
737  if (prob.colstat_)
738    presolve_check_nbasic(&prob) ;
739 
740//check_djs(&prob);
741  presolve_check_reduced_costs(&prob) ;
742#endif
743 
744 
745  while (paction) {
746#if     PRESOLVE_DEBUG
747    printf("POSTSOLVING %s\n", paction->name());
748#endif
749
750    paction->postsolve(&prob);
751   
752#if     PRESOLVE_DEBUG
753    if (prob.colstat_)
754      presolve_check_nbasic(&prob) ;
755#endif
756    paction = paction->next;
757#if     PRESOLVE_DEBUG
758//  check_djs(&prob);
759    presolve_check_reduced_costs(&prob) ;
760#endif
761  }   
762 
763#if     0 && PRESOLVE_DEBUG
764  for (i=0; i<ncols0; i++) {
765    if (!cdone[i]) {
766      printf("!cdone[%d]\n", i);
767      abort();
768    }
769  }
770 
771  for (i=0; i<nrows0; i++) {
772    if (!rdone[i]) {
773      printf("!rdone[%d]\n", i);
774      abort();
775    }
776  }
777 
778 
779  for (i=0; i<ncols0; i++) {
780    if (sol[i] < -1e10 || sol[i] > 1e10)
781      printf("!!!%d %g\n", i, sol[i]);
782   
783  }
784 
785 
786#endif
787 
788#if     0 && PRESOLVE_DEBUG
789  // debug check:  make sure we ended up with same original matrix
790  {
791    int identical = 1;
792   
793    for (int i=0; i<ncols0; i++) {
794      PRESOLVEASSERT(hincol[i] == &prob->mcstrt0[i+1] - &prob->mcstrt0[i]);
795      CoinBigIndex kcs0 = &prob->mcstrt0[i];
796      CoinBigIndex kcs = mcstrt[i];
797      int n = hincol[i];
798      for (int k=0; k<n; k++) {
799        CoinBigIndex k1 = presolve_find_row1(&prob->hrow0[kcs0+k], kcs, kcs+n, hrow);
800
801        if (k1 == kcs+n) {
802          printf("ROW %d NOT IN COL %d\n", &prob->hrow0[kcs0+k], i);
803          abort();
804        }
805
806        if (colels[k1] != &prob->dels0[kcs0+k])
807          printf("BAD COLEL[%d %d %d]:  %g\n",
808                 k1, i, &prob->hrow0[kcs0+k], colels[k1] - &prob->dels0[kcs0+k]);
809
810        if (kcs0+k != k1)
811          identical=0;
812      }
813    }
814    printf("identical? %d\n", identical);
815  }
816#endif
817}
818
819
820
821
822
823
824
825
826static inline double getTolerance(const ClpSimplex  *si, ClpDblParam key)
827{
828  double tol;
829  if (! si->getDblParam(key, tol)) {
830    CoinPresolveAction::throwCoinError("getDblParam failed",
831                                      "CoinPrePostsolveMatrix::CoinPrePostsolveMatrix");
832  }
833  return (tol);
834}
835
836
837// Assumptions:
838// 1. nrows>=m.getNumRows()
839// 2. ncols>=m.getNumCols()
840//
841// In presolve, these values are equal.
842// In postsolve, they may be inequal, since the reduced problem
843// may be smaller, but we need room for the large problem.
844// ncols may be larger than si.getNumCols() in postsolve,
845// this at that point si will be the reduced problem,
846// but we need to reserve enough space for the original problem.
847CoinPrePostsolveMatrix::CoinPrePostsolveMatrix(const ClpSimplex * si,
848                                             int ncols_in,
849                                             int nrows_in,
850                                             CoinBigIndex nelems_in) 
851  : ncols_(si->getNumCols()),
852    nrows_(si->getNumRows()),
853    nelems_(si->getNumElements()),
854    ncols0_(ncols_in),
855    nrows0_(nrows_in),
856    mcstrt_(new CoinBigIndex[ncols_in+1]),
857    hincol_(new int[ncols_in+1]),
858    hrow_  (new int   [2*nelems_in]),
859    colels_(new double[2*nelems_in]),
860    cost_(new double[ncols_in]),
861    clo_(new double[ncols_in]),
862    cup_(new double[ncols_in]),
863    rlo_(new double[nrows_in]),
864    rup_(new double[nrows_in]),
865    originalColumn_(new int[ncols_in]),
866    originalRow_(new int[nrows_in]),
867    ztolzb_(getTolerance(si, ClpPrimalTolerance)),
868    ztoldj_(getTolerance(si, ClpDualTolerance)),
869    maxmin_(si->getObjSense()),
870    sol_(NULL),
871    rowduals_(NULL),
872    acts_(NULL),
873    rcosts_(NULL),
874    colstat_(NULL),
875    rowstat_(NULL),
876    handler_(NULL),
877    defaultHandler_(false),
878    messages_()
879
880{
881  si->getDblParam(ClpObjOffset,originalOffset_);
882  int ncols = si->getNumCols();
883  int nrows = si->getNumRows();
884
885  setMessageHandler(si->messageHandler()) ;
886
887  ClpDisjointCopyN(si->getColLower(), ncols, clo_);
888  ClpDisjointCopyN(si->getColUpper(), ncols, cup_);
889  //ClpDisjointCopyN(si->getObjCoefficients(), ncols, cost_);
890  double offset;
891  ClpDisjointCopyN(si->objectiveAsObject()->gradient(si,si->getColSolution(),offset,true), ncols, cost_);
892  ClpDisjointCopyN(si->getRowLower(), nrows,  rlo_);
893  ClpDisjointCopyN(si->getRowUpper(), nrows,  rup_);
894  int i;
895  for (i=0;i<ncols_in;i++) 
896    originalColumn_[i]=i;
897  for (i=0;i<nrows_in;i++) 
898    originalRow_[i]=i;
899  sol_=NULL;
900  rowduals_=NULL;
901  acts_=NULL;
902
903  rcosts_=NULL;
904  colstat_=NULL;
905  rowstat_=NULL;
906}
907
908// I am not familiar enough with CoinPackedMatrix to be confident
909// that I will implement a row-ordered version of toColumnOrderedGapFree
910// properly.
911static bool isGapFree(const CoinPackedMatrix& matrix)
912{
913  const CoinBigIndex * start = matrix.getVectorStarts();
914  const int * length = matrix.getVectorLengths();
915  int i;
916  for (i = matrix.getSizeVectorLengths() - 1; i >= 0; --i) {
917    if (start[i+1] - start[i] != length[i])
918      break;
919  }
920  return (! (i >= 0));
921}
922#if     PRESOLVE_DEBUG
923static void matrix_bounds_ok(const double *lo, const double *up, int n)
924{
925  int i;
926  for (i=0; i<n; i++) {
927    PRESOLVEASSERT(lo[i] <= up[i]);
928    PRESOLVEASSERT(lo[i] < PRESOLVE_INF);
929    PRESOLVEASSERT(-PRESOLVE_INF < up[i]);
930  }
931}
932#endif
933CoinPresolveMatrix::CoinPresolveMatrix(int ncols0_in,
934                                     double maxmin,
935                                     // end prepost members
936
937                                     ClpSimplex * si,
938
939                                     // rowrep
940                                     int nrows_in,
941                                     CoinBigIndex nelems_in,
942                               bool doStatus,
943                               double nonLinearValue) :
944
945  CoinPrePostsolveMatrix(si,
946                        ncols0_in, nrows_in, nelems_in),
947  clink_(new presolvehlink[ncols0_in+1]),
948  rlink_(new presolvehlink[nrows_in+1]),
949
950  dobias_(0.0),
951
952
953  // temporary init
954  mrstrt_(new CoinBigIndex[nrows_in+1]),
955  hinrow_(new int[nrows_in+1]),
956  rowels_(new double[2*nelems_in]),
957  hcol_(new int[2*nelems_in]),
958  integerType_(new unsigned char[ncols0_in]),
959  tuning_(false),
960  startTime_(0.0),
961  feasibilityTolerance_(0.0),
962  status_(-1),
963  colsToDo_(new int [ncols0_in]),
964  numberColsToDo_(0),
965  nextColsToDo_(new int[ncols0_in]),
966  numberNextColsToDo_(0),
967  rowsToDo_(new int [nrows_in]),
968  numberRowsToDo_(0),
969  nextRowsToDo_(new int[nrows_in]),
970  numberNextRowsToDo_(0)
971
972{
973  const int bufsize = 2*nelems_in;
974
975  nrows_ = si->getNumRows() ;
976
977  // Set up change bits etc
978  rowChanged_ = new unsigned char[nrows_];
979  memset(rowChanged_,0,nrows_);
980  colChanged_ = new unsigned char[ncols_];
981  memset(colChanged_,0,ncols_);
982  CoinPackedMatrix * m = si->matrix();
983
984  // The coefficient matrix is a big hunk of stuff.
985  // Do the copy here to try to avoid running out of memory.
986
987  const CoinBigIndex * start = m->getVectorStarts();
988  const int * length = m->getVectorLengths();
989  const int * row = m->getIndices();
990  const double * element = m->getElements();
991  int icol,nel=0;
992  mcstrt_[0]=0;
993  for (icol=0;icol<ncols_;icol++) {
994    int j;
995    for (j=start[icol];j<start[icol]+length[icol];j++) {
996      hrow_[nel]=row[j];
997      colels_[nel++]=element[j];
998    }
999    mcstrt_[icol+1]=nel;
1000  }
1001  assert(mcstrt_[ncols_] == nelems_);
1002  ClpDisjointCopyN(m->getVectorLengths(),ncols_,  hincol_);
1003
1004  // same thing for row rep
1005  m = new CoinPackedMatrix();
1006  m->reverseOrderedCopyOf(*si->matrix());
1007  m->removeGaps();
1008
1009
1010  ClpDisjointCopyN(m->getVectorStarts(),  nrows_,  mrstrt_);
1011  mrstrt_[nrows_] = nelems_;
1012  ClpDisjointCopyN(m->getVectorLengths(), nrows_,  hinrow_);
1013  ClpDisjointCopyN(m->getIndices(),       nelems_, hcol_);
1014  ClpDisjointCopyN(m->getElements(),      nelems_, rowels_);
1015
1016  delete m;
1017  if (si->integerInformation()) {
1018    memcpy(integerType_,si->integerInformation(),ncols_*sizeof(char));
1019  } else {
1020    ClpFillN<unsigned char>(integerType_, ncols_, (unsigned char) 0);
1021  }
1022
1023#ifndef NO_RTTI
1024    ClpQuadraticObjective * quadraticObj = (dynamic_cast< ClpQuadraticObjective*>(si->objectiveAsObject()));
1025#else
1026    ClpQuadraticObjective * quadraticObj = NULL;
1027    if (si->objectiveAsObject()->type()==2)
1028      quadraticObj = (static_cast< ClpQuadraticObjective*>(si->objectiveAsObject()));
1029#endif
1030  // Set up prohibited bits if needed
1031  if (nonLinearValue) {
1032    anyProhibited_ = true;
1033    for (icol=0;icol<ncols_;icol++) {
1034      int j;
1035      bool nonLinearColumn = false;
1036      if (cost_[icol]==nonLinearValue)
1037        nonLinearColumn=true;
1038      for (j=mcstrt_[icol];j<mcstrt_[icol+1];j++) {
1039        if (colels_[j]==nonLinearValue) {
1040          nonLinearColumn=true;
1041          setRowProhibited(hrow_[j]);
1042        }
1043      }
1044      if (nonLinearColumn)
1045        setColProhibited(icol);
1046    }
1047  } else if (quadraticObj) {
1048    CoinPackedMatrix * quadratic = quadraticObj->quadraticObjective();
1049    //const int * columnQuadratic = quadratic->getIndices();
1050    //const CoinBigIndex * columnQuadraticStart = quadratic->getVectorStarts();
1051    const int * columnQuadraticLength = quadratic->getVectorLengths();
1052    //double * quadraticElement = quadratic->getMutableElements();
1053    int numberColumns = quadratic->getNumCols();
1054    anyProhibited_ = true;
1055    for (int iColumn=0;iColumn<numberColumns;iColumn++) {
1056      if (columnQuadraticLength[iColumn]) {
1057        setColProhibited(iColumn);
1058        //printf("%d prohib\n",iColumn);
1059      }
1060    }
1061  } else {
1062    anyProhibited_ = false;
1063  }
1064
1065  if (doStatus) {
1066    // allow for status and solution
1067    sol_ = new double[ncols_];
1068    memcpy(sol_,si->primalColumnSolution(),ncols_*sizeof(double));;
1069    acts_ = new double [nrows_];
1070    memcpy(acts_,si->primalRowSolution(),nrows_*sizeof(double));
1071    if (!si->statusArray())
1072      si->createStatus();
1073    colstat_ = new unsigned char [nrows_+ncols_];
1074    memcpy(colstat_,si->statusArray(),
1075           (nrows_+ncols_)*sizeof(unsigned char));
1076    rowstat_ = colstat_+ncols_;
1077  }
1078
1079  // the original model's fields are now unneeded - free them
1080 
1081  si->resize(0,0);
1082
1083#if     PRESOLVE_DEBUG
1084  matrix_bounds_ok(rlo_, rup_, nrows_);
1085  matrix_bounds_ok(clo_, cup_, ncols_);
1086#endif
1087
1088#if 0
1089  for (i=0; i<nrows; ++i)
1090    printf("NR: %6d\n", hinrow[i]);
1091  for (int i=0; i<ncols; ++i)
1092    printf("NC: %6d\n", hincol[i]);
1093#endif
1094
1095  presolve_make_memlists(mcstrt_, hincol_, clink_, ncols_);
1096  presolve_make_memlists(mrstrt_, hinrow_, rlink_, nrows_);
1097
1098  // this allows last col/row to expand up to bufsize-1 (22);
1099  // this must come after the calls to presolve_prefix
1100  mcstrt_[ncols_] = bufsize-1;
1101  mrstrt_[nrows_] = bufsize-1;
1102
1103#if     PRESOLVE_CONSISTENCY
1104//consistent(false);
1105  presolve_consistent(this,false) ;
1106#endif
1107}
1108
1109void CoinPresolveMatrix::update_model(ClpSimplex * si,
1110                                     int nrows0,
1111                                     int ncols0,
1112                                     CoinBigIndex nelems0)
1113{
1114  si->loadProblem(ncols_, nrows_, mcstrt_, hrow_, colels_, hincol_,
1115                 clo_, cup_, cost_, rlo_, rup_);
1116  delete [] si->integerInformation();
1117  int numberIntegers=0;
1118  for (int i=0; i<ncols_; i++) {
1119    if (integerType_[i])
1120      numberIntegers++;
1121  }
1122  if (numberIntegers) 
1123    si->copyInIntegerInformation((const char *) integerType_);
1124  else
1125    si->copyInIntegerInformation(NULL);
1126
1127#if     PRESOLVE_SUMMARY
1128  printf("NEW NCOL/NROW/NELS:  %d(-%d) %d(-%d) %d(-%d)\n",
1129         ncols_, ncols0-ncols_,
1130         nrows_, nrows0-nrows_,
1131         si->getNumElements(), nelems0-si->getNumElements());
1132#endif
1133  si->setDblParam(ClpObjOffset,originalOffset_-dobias_);
1134
1135}
1136
1137
1138
1139
1140
1141
1142
1143
1144
1145
1146
1147////////////////  POSTSOLVE
1148
1149CoinPostsolveMatrix::CoinPostsolveMatrix(ClpSimplex*  si,
1150                                       int ncols0_in,
1151                                       int nrows0_in,
1152                                       CoinBigIndex nelems0,
1153                                   
1154                                       double maxmin,
1155                                       // end prepost members
1156
1157                                       double *sol_in,
1158                                       double *acts_in,
1159
1160                                       unsigned char *colstat_in,
1161                                       unsigned char *rowstat_in) :
1162  CoinPrePostsolveMatrix(si,
1163                        ncols0_in, nrows0_in, nelems0),
1164
1165  free_list_(0),
1166  // link, free_list, maxlink
1167  maxlink_(2*nelems0),
1168  link_(new int[/*maxlink*/ 2*nelems0]),
1169     
1170  cdone_(new char[ncols0_]),
1171  rdone_(new char[nrows0_in])
1172
1173{
1174  bulk0_ = maxlink_ ;
1175  nrows_ = si->getNumRows() ;
1176  ncols_ = si->getNumCols() ;
1177
1178  sol_=sol_in;
1179  rowduals_=NULL;
1180  acts_=acts_in;
1181
1182  rcosts_=NULL;
1183  colstat_=colstat_in;
1184  rowstat_=rowstat_in;
1185
1186  // this is the *reduced* model, which is probably smaller
1187  int ncols1 = ncols_ ;
1188  int nrows1 = nrows_ ;
1189
1190  const CoinPackedMatrix * m = si->matrix();
1191
1192  if (! isGapFree(*m)) {
1193    CoinPresolveAction::throwCoinError("Matrix not gap free",
1194                                      "CoinPostsolveMatrix");
1195  }
1196
1197  const CoinBigIndex nelemsr = m->getNumElements();
1198
1199  ClpDisjointCopyN(m->getVectorStarts(), ncols1, mcstrt_);
1200  CoinZeroN(mcstrt_+ncols1,ncols0_-ncols1);
1201  mcstrt_[ncols1] = nelems0;    // ??    (should point to end of bulk store   -- lh --)
1202  ClpDisjointCopyN(m->getVectorLengths(),ncols1,  hincol_);
1203  ClpDisjointCopyN(m->getIndices(),      nelemsr, hrow_);
1204  ClpDisjointCopyN(m->getElements(),     nelemsr, colels_);
1205
1206
1207#if     0 && PRESOLVE_DEBUG
1208  presolve_check_costs(model, &colcopy);
1209#endif
1210
1211  // This determines the size of the data structure that contains
1212  // the matrix being postsolved.  Links are taken from the free_list
1213  // to recreate matrix entries that were presolved away,
1214  // and links are added to the free_list when entries created during
1215  // presolve are discarded.  There is never a need to gc this list.
1216  // Naturally, it should contain
1217  // exactly nelems0 entries "in use" when postsolving is done,
1218  // but I don't know whether the matrix could temporarily get
1219  // larger during postsolving.  Substitution into more than two
1220  // rows could do that, in principle.  I am being very conservative
1221  // here by reserving much more than the amount of space I probably need.
1222  // If this guess is wrong, check_free_list may be called.
1223  //  int bufsize = 2*nelems0;
1224
1225  memset(cdone_, -1, ncols0_);
1226  memset(rdone_, -1, nrows0_);
1227
1228  rowduals_ = new double[nrows0_];
1229  ClpDisjointCopyN(si->getRowPrice(), nrows1, rowduals_);
1230
1231  rcosts_ = new double[ncols0_];
1232  ClpDisjointCopyN(si->getReducedCost(), ncols1, rcosts_);
1233  if (maxmin<0.0) {
1234    // change so will look as if minimize
1235    int i;
1236    for (i=0;i<nrows1;i++)
1237      rowduals_[i] = - rowduals_[i];
1238    for (i=0;i<ncols1;i++) {
1239      rcosts_[i] = - rcosts_[i];
1240    }
1241  }
1242
1243  //ClpDisjointCopyN(si->getRowUpper(), nrows1, rup_);
1244  //ClpDisjointCopyN(si->getRowLower(), nrows1, rlo_);
1245
1246  ClpDisjointCopyN(si->getColSolution(), ncols1, sol_);
1247  si->setDblParam(ClpObjOffset,originalOffset_);
1248
1249  for (int j=0; j<ncols1; j++) {
1250    CoinBigIndex kcs = mcstrt_[j];
1251    CoinBigIndex kce = kcs + hincol_[j];
1252    for (CoinBigIndex k=kcs; k<kce; ++k) {
1253      link_[k] = k+1;
1254    }
1255    link_[kce-1] = NO_LINK ;
1256  }
1257  {
1258    int ml = maxlink_;
1259    for (CoinBigIndex k=nelemsr; k<ml; ++k)
1260      link_[k] = k+1;
1261    if (ml)
1262      link_[ml-1] = NO_LINK;
1263  }
1264  free_list_ = nelemsr;
1265# if PRESOLVE_DEBUG || PRESOLVE_CONSISTENCY
1266/*
1267  These are used to track the action of postsolve transforms during debugging.
1268*/
1269  CoinFillN(cdone_,ncols1,PRESENT_IN_REDUCED) ;
1270  CoinZeroN(cdone_+ncols1,ncols0_in-ncols1) ;
1271  CoinFillN(rdone_,nrows1,PRESENT_IN_REDUCED) ;
1272  CoinZeroN(rdone_+nrows1,nrows0_in-nrows1) ;
1273# endif
1274}
1275/* This is main part of Presolve */
1276ClpSimplex * 
1277ClpPresolve::gutsOfPresolvedModel(ClpSimplex * originalModel,
1278                                  double feasibilityTolerance,
1279                                  bool keepIntegers,
1280                                  int numberPasses,
1281                                  bool dropNames)
1282{
1283  ncols_ = originalModel->getNumCols();
1284  nrows_ = originalModel->getNumRows();
1285  nelems_ = originalModel->getNumElements();
1286  numberPasses_ = numberPasses;
1287
1288  double maxmin = originalModel->getObjSense();
1289  originalModel_ = originalModel;
1290  delete [] originalColumn_;
1291  originalColumn_ = new int[ncols_];
1292  delete [] originalRow_;
1293  originalRow_ = new int[nrows_];
1294
1295  // result is 0 - okay, 1 infeasible, -1 go round again, 2 - original model
1296  int result = -1;
1297 
1298  // User may have deleted - its their responsibility
1299  presolvedModel_=NULL;
1300  // Messages
1301  CoinMessages messages = CoinMessage(originalModel->messages().language());
1302  while (result==-1) {
1303
1304    // make new copy
1305    if (saveFile_=="") {
1306      delete presolvedModel_;
1307      presolvedModel_ = new ClpSimplex(*originalModel);
1308    } else {
1309      presolvedModel_=originalModel;
1310    }
1311    presolvedModel_->dropNames();
1312
1313    // drop integer information if wanted
1314    if (!keepIntegers)
1315      presolvedModel_->deleteIntegerInformation();
1316
1317   
1318    CoinPresolveMatrix prob(ncols_,
1319                        maxmin,
1320                        presolvedModel_,
1321                        nrows_, nelems_,true,nonLinearValue_);
1322    // See if we want statistics
1323    if ((presolveActions_&0x80000000)!=0)
1324      prob.statistics();
1325    // make sure row solution correct
1326    {
1327      double *colels    = prob.colels_;
1328      int *hrow         = prob.hrow_;
1329      CoinBigIndex *mcstrt              = prob.mcstrt_;
1330      int *hincol               = prob.hincol_;
1331      int ncols         = prob.ncols_;
1332     
1333     
1334      double * csol = prob.sol_;
1335      double * acts = prob.acts_;
1336      int nrows = prob.nrows_;
1337
1338      int colx;
1339
1340      memset(acts,0,nrows*sizeof(double));
1341     
1342      for (colx = 0; colx < ncols; ++colx) {
1343        double solutionValue = csol[colx];
1344        for (int i=mcstrt[colx]; i<mcstrt[colx]+hincol[colx]; ++i) {
1345          int row = hrow[i];
1346          double coeff = colels[i];
1347          acts[row] += solutionValue*coeff;
1348        }
1349      }
1350    }
1351
1352    // move across feasibility tolerance
1353    prob.feasibilityTolerance_ = feasibilityTolerance;
1354
1355    // Do presolve
1356    paction_ = presolve(&prob);
1357
1358    result =0; 
1359
1360    if (prob.status_==0&&paction_) {
1361      // Looks feasible but double check to see if anything slipped through
1362      int n             = prob.ncols_;
1363      double * lo = prob.clo_;
1364      double * up = prob.cup_;
1365      int i;
1366     
1367      for (i=0;i<n;i++) {
1368        if (up[i]<lo[i]) {
1369          if (up[i]<lo[i]-1.0e-8) {
1370            // infeasible
1371            prob.status_=1;
1372          } else {
1373            up[i]=lo[i];
1374          }
1375        }
1376      }
1377     
1378      n = prob.nrows_;
1379      lo = prob.rlo_;
1380      up = prob.rup_;
1381
1382      for (i=0;i<n;i++) {
1383        if (up[i]<lo[i]) {
1384          if (up[i]<lo[i]-1.0e-8) {
1385            // infeasible
1386            prob.status_=1;
1387          } else {
1388            up[i]=lo[i];
1389          }
1390        }
1391      }
1392    }
1393    if (prob.status_==0&&paction_) {
1394      // feasible
1395   
1396      prob.update_model(presolvedModel_, nrows_, ncols_, nelems_);
1397      // copy status and solution
1398      memcpy(presolvedModel_->primalColumnSolution(),
1399             prob.sol_,prob.ncols_*sizeof(double));
1400      memcpy(presolvedModel_->primalRowSolution(),
1401             prob.acts_,prob.nrows_*sizeof(double));
1402      memcpy(presolvedModel_->statusArray(),
1403             prob.colstat_,prob.ncols_*sizeof(unsigned char));
1404      memcpy(presolvedModel_->statusArray()+prob.ncols_,
1405             prob.rowstat_,prob.nrows_*sizeof(unsigned char));
1406      delete [] prob.sol_;
1407      delete [] prob.acts_;
1408      delete [] prob.colstat_;
1409      prob.sol_=NULL;
1410      prob.acts_=NULL;
1411      prob.colstat_=NULL;
1412     
1413      int ncolsNow = presolvedModel_->getNumCols();
1414      memcpy(originalColumn_,prob.originalColumn_,ncolsNow*sizeof(int));
1415#ifndef NO_RTTI
1416      ClpQuadraticObjective * quadraticObj = (dynamic_cast< ClpQuadraticObjective*>(originalModel->objectiveAsObject()));
1417#else
1418      ClpQuadraticObjective * quadraticObj = NULL;
1419      if (originalModel->objectiveAsObject()->type()==2)
1420        quadraticObj = (static_cast< ClpQuadraticObjective*>(originalModel->objectiveAsObject()));
1421#endif
1422      if (quadraticObj) {
1423        // set up for subset
1424        char * mark = new char [ncols_];
1425        memset(mark,0,ncols_);
1426        CoinPackedMatrix * quadratic = quadraticObj->quadraticObjective();
1427        //const int * columnQuadratic = quadratic->getIndices();
1428        //const CoinBigIndex * columnQuadraticStart = quadratic->getVectorStarts();
1429        const int * columnQuadraticLength = quadratic->getVectorLengths();
1430        //double * quadraticElement = quadratic->getMutableElements();
1431        int numberColumns = quadratic->getNumCols();
1432        ClpQuadraticObjective * newObj = new ClpQuadraticObjective(*quadraticObj,
1433                                                                   ncolsNow,
1434                                                                   originalColumn_);
1435        // and modify linear and check
1436        double * linear = newObj->linearObjective();
1437        memcpy(linear,presolvedModel_->objective(),ncolsNow*sizeof(double));
1438        int iColumn;
1439        for ( iColumn=0;iColumn<numberColumns;iColumn++) {
1440          if (columnQuadraticLength[iColumn])
1441            mark[iColumn]=1;
1442        }
1443        // and new
1444        quadratic = newObj->quadraticObjective();
1445        columnQuadraticLength = quadratic->getVectorLengths();
1446        int numberColumns2 = quadratic->getNumCols();
1447        for ( iColumn=0;iColumn<numberColumns2;iColumn++) {
1448          if (columnQuadraticLength[iColumn])
1449            mark[originalColumn_[iColumn]]=0;
1450        }
1451        presolvedModel_->setObjective(newObj);
1452        delete newObj;
1453        // final check
1454        for ( iColumn=0;iColumn<numberColumns;iColumn++) 
1455          if (mark[iColumn])
1456            printf("Quadratic column %d modified - may be okay\n",iColumn);
1457        delete [] mark;
1458      }
1459      delete [] prob.originalColumn_;
1460      prob.originalColumn_=NULL;
1461      int nrowsNow = presolvedModel_->getNumRows();
1462      memcpy(originalRow_,prob.originalRow_,nrowsNow*sizeof(int));
1463      delete [] prob.originalRow_;
1464      prob.originalRow_=NULL;
1465      if (!dropNames&&originalModel->lengthNames()) {
1466        // Redo names
1467        int iRow;
1468        std::vector<std::string> rowNames;
1469        rowNames.reserve(nrowsNow);
1470        for (iRow=0;iRow<nrowsNow;iRow++) {
1471          int kRow = originalRow_[iRow];
1472          rowNames.push_back(originalModel->rowName(kRow));
1473        }
1474     
1475        int iColumn;
1476        std::vector<std::string> columnNames;
1477        columnNames.reserve(ncolsNow);
1478        for (iColumn=0;iColumn<ncolsNow;iColumn++) {
1479          int kColumn = originalColumn_[iColumn];
1480          columnNames.push_back(originalModel->columnName(kColumn));
1481        }
1482        presolvedModel_->copyNames(rowNames,columnNames);
1483      }
1484      // now clean up integer variables.  This can modify original
1485      int i;
1486      const char * information = presolvedModel_->integerInformation();
1487      if (information) {
1488        int numberChanges=0;
1489        double * lower0 = originalModel_->columnLower();
1490        double * upper0 = originalModel_->columnUpper();
1491        double * lower = presolvedModel_->columnLower();
1492        double * upper = presolvedModel_->columnUpper();
1493        for (i=0;i<ncolsNow;i++) {
1494          if (!information[i])
1495            continue;
1496          int iOriginal = originalColumn_[i];
1497          double lowerValue0 = lower0[iOriginal];
1498          double upperValue0 = upper0[iOriginal];
1499          double lowerValue = ceil(lower[i]-1.0e-5);
1500          double upperValue = floor(upper[i]+1.0e-5);
1501          lower[i]=lowerValue;
1502          upper[i]=upperValue;
1503          if (lowerValue>upperValue) {
1504            numberChanges++;
1505            presolvedModel_->messageHandler()->message(COIN_PRESOLVE_COLINFEAS,
1506                                                       messages)
1507                                                         <<iOriginal
1508                                                         <<lowerValue
1509                                                         <<upperValue
1510                                                         <<CoinMessageEol;
1511            result=1;
1512          } else {
1513            if (lowerValue>lowerValue0+1.0e-8) {
1514              lower0[iOriginal] = lowerValue;
1515              numberChanges++;
1516            }
1517            if (upperValue<upperValue0-1.0e-8) {
1518              upper0[iOriginal] = upperValue;
1519              numberChanges++;
1520            }
1521          }       
1522        }
1523        if (numberChanges) {
1524          presolvedModel_->messageHandler()->message(COIN_PRESOLVE_INTEGERMODS,
1525                                                     messages)
1526                                                       <<numberChanges
1527                                                       <<CoinMessageEol;
1528          if (!result) {
1529            result = -1; // round again
1530          }
1531        }
1532      }
1533    } else if (prob.status_) {
1534      // infeasible or unbounded
1535      result=1;
1536    } else {
1537      // no changes - model needs restoring after Lou's changes
1538      if (saveFile_=="") {
1539        delete presolvedModel_;
1540        presolvedModel_ = new ClpSimplex(*originalModel);
1541      } else {
1542        presolvedModel_=originalModel;
1543      }
1544      presolvedModel_->dropNames();
1545     
1546      // drop integer information if wanted
1547      if (!keepIntegers)
1548        presolvedModel_->deleteIntegerInformation();
1549      result=2;
1550    }
1551  }
1552  if (result==0||result==2) {
1553    int nrowsAfter = presolvedModel_->getNumRows();
1554    int ncolsAfter = presolvedModel_->getNumCols();
1555    CoinBigIndex nelsAfter = presolvedModel_->getNumElements();
1556    presolvedModel_->messageHandler()->message(COIN_PRESOLVE_STATS,
1557                                               messages)
1558                                                 <<nrowsAfter<< -(nrows_ - nrowsAfter)
1559                                                 << ncolsAfter<< -(ncols_ - ncolsAfter)
1560                                                 <<nelsAfter<< -(nelems_ - nelsAfter)
1561                                                 <<CoinMessageEol;
1562  } else {
1563    destroyPresolve();
1564    delete presolvedModel_;
1565    presolvedModel_=NULL;
1566  }
1567  return presolvedModel_;
1568}
1569
1570
Note: See TracBrowser for help on using the repository browser.