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

Last change on this file since 1551 was 1551, checked in by mjs, 9 years ago

Rerun formatting after JJF patches (1533 and 1535).

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