source: trunk/ClpPresolve.cpp @ 446

Last change on this file since 446 was 446, checked in by lou-oss, 15 years ago

Compatibility w/ CoinPresolve? modifications.

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