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

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

allow for network bobble and lightweight presolve

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