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

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

fix bug when presolved model infeasible but original feasible.
also allow names in presolve to file.

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