source: trunk/ClpPresolve.cpp @ 458

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

So can switch off dual part of presolve

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