source: trunk/ClpPresolve.cpp @ 393

Last change on this file since 393 was 393, checked in by forrest, 16 years ago

Some quadratic stuff

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