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

Last change on this file since 1660 was 1660, checked in by forrest, 9 years ago

for user hooks and ranging/parametrics

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