source: trunk/ClpPresolve.cpp @ 582

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

going for robustness

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