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

Last change on this file since 1665 was 1665, checked in by lou, 9 years ago

Add EPL license notice in src.

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