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

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

allow to distinguish between infeasible and unbounded in presolve

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