source: trunk/ClpPresolve.cpp @ 450

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

For Lou's changes

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