source: trunk/ClpPresolve.cpp @ 564

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

marginally faster if not CoinPackedMatrix?

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