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

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

add debug info and take out possible memory leak

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