source: trunk/Clp/src/ClpPresolve.cpp @ 1834

Last change on this file since 1834 was 1834, checked in by lou, 8 years ago

Disable some asserts enabled by PRESOLVE_DEBUG. Clp appears to return
incorrect status at start of postsolve, hence checks fail.

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