source: trunk/Clp/examples/sprint2.cpp

Last change on this file was 1937, checked in by forrest, 5 years ago

fixes so sprint will work on unsuitable problems and testGub have no memory leaks

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 6.5 KB
Line 
1/* $Id: sprint2.cpp 1937 2013-04-09 10:32:15Z forrest $ */
2// Copyright (C) 2003, International Business Machines
3// Corporation and others.  All Rights Reserved.
4// This code is licensed under the terms of the Eclipse Public License (EPL).
5
6#include "ClpSimplex.hpp"
7#include "ClpPresolve.hpp"
8#include "CoinSort.hpp"
9#include <iomanip>
10
11int main(int argc, const char *argv[])
12{
13     ClpSimplex  model;
14     int status;
15     // Keep names
16     if (argc < 2) {
17          status = model.readMps("small.mps", true);
18     } else {
19          status = model.readMps(argv[1], true);
20     }
21     if (status)
22          exit(10);
23     /*
24       This driver implements the presolve variation of Sprint.
25       This assumes we can get feasible easily
26     */
27
28     int numberRows = model.numberRows();
29     int numberColumns = model.numberColumns();
30
31     // We will need arrays to choose variables.  These are too big but ..
32     double * weight = new double [numberRows+numberColumns];
33     int * sort = new int [numberRows+numberColumns];
34
35     double * columnLower = model.columnLower();
36     double * columnUpper = model.columnUpper();
37     double * saveLower = new double [numberColumns];
38     memcpy(saveLower, columnLower, numberColumns * sizeof(double));
39     double * saveUpper = new double [numberColumns];
40     memcpy(saveUpper, columnUpper, numberColumns * sizeof(double));
41     double * solution = model.primalColumnSolution();
42     // Fix in some magical way so remaining problem is easy
43#if 0
44     // This is from a real-world problem
45     for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
46          char firstCharacter = model.columnName(iColumn)[0];
47          if (firstCharacter == 'F' || firstCharacter == 'P'
48                    || firstCharacter == 'L' || firstCharacter == 'T') {
49               columnUpper[iColumn] = columnLower[iColumn];
50          }
51     }
52#else
53     double * obj = model.objective();
54     double * saveObj = new double [numberColumns];
55     memcpy(saveObj, obj, numberColumns * sizeof(double));
56     memset(obj, 0, numberColumns * sizeof(double));
57     model.dual();
58     memcpy(obj, saveObj, numberColumns * sizeof(double));
59     delete [] saveObj;
60     for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
61         if (solution[iColumn]<columnLower[iColumn]+1.0e-8) {
62             columnUpper[iColumn] = columnLower[iColumn];
63       }
64     }
65#endif
66
67     // Just do this number of passes
68     int maxPass = 100;
69     int iPass;
70     double lastObjective = 1.0e31;
71
72     // Just take this number of columns in small problem
73     int smallNumberColumns = 3 * numberRows;
74     // To stop seg faults on unsuitable problems
75     smallNumberColumns = CoinMin(smallNumberColumns,numberColumns);
76     // And we want number of rows to be this
77     int smallNumberRows = numberRows / 4;
78
79     for (iPass = 0; iPass < maxPass; iPass++) {
80          printf("Start of pass %d\n", iPass);
81          ClpSimplex * model2;
82          ClpPresolve pinfo;
83          int numberPasses = 1; // can change this
84          model2 = pinfo.presolvedModel(model, 1.0e-8, false, numberPasses, false);
85          if (!model2) {
86               fprintf(stdout, "ClpPresolve says %s is infeasible with tolerance of %g\n",
87                       argv[1], 1.0e-8);
88               // model was infeasible - maybe try again with looser tolerances
89               model2 = pinfo.presolvedModel(model, 1.0e-7, false, numberPasses, false);
90               if (!model2) {
91                    fprintf(stdout, "ClpPresolve says %s is infeasible with tolerance of %g\n",
92                            argv[1], 1.0e-7);
93                    exit(2);
94               }
95          }
96          // change factorization frequency from 200
97          model2->setFactorizationFrequency(100 + model2->numberRows() / 50);
98          model2->primal();
99          pinfo.postsolve(true);
100
101          // adjust smallNumberColumns if necessary
102          if (iPass) {
103               double ratio = ((double) smallNumberRows) / ((double) model2->numberRows());
104               smallNumberColumns = (int)(smallNumberColumns * ratio);
105               // deal with pathological case
106               smallNumberColumns = CoinMax(smallNumberColumns,0);
107          }
108          delete model2;
109          /* After this postsolve model should be optimal.
110             We can use checkSolution and test feasibility */
111          model.checkSolution();
112          if (model.numberDualInfeasibilities() ||
113                    model.numberPrimalInfeasibilities())
114               printf("%g dual %g(%d) Primal %g(%d)\n",
115                      model.objectiveValue(),
116                      model.sumDualInfeasibilities(),
117                      model.numberDualInfeasibilities(),
118                      model.sumPrimalInfeasibilities(),
119                      model.numberPrimalInfeasibilities());
120          // Put back true bounds
121          memcpy(columnLower, saveLower, numberColumns * sizeof(double));
122          memcpy(columnUpper, saveUpper, numberColumns * sizeof(double));
123          if ((model.objectiveValue() > lastObjective - 1.0e-7 && iPass > 5) ||
124                    iPass == maxPass - 1) {
125               break; // finished
126          } else {
127               lastObjective = model.objectiveValue();
128               // now massage weight so all basic in plus good djs
129               const double * djs = model.dualColumnSolution();
130               for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
131                    double dj = djs[iColumn];
132                    double value = solution[iColumn];
133                    if (model.getStatus(iColumn) == ClpSimplex::basic)
134                         dj = -1.0e50;
135                    else if (dj < 0.0 && value < columnUpper[iColumn])
136                         dj = dj;
137                    else if (dj > 0.0 && value > columnLower[iColumn])
138                         dj = -dj;
139                    else if (columnUpper[iColumn] > columnLower[iColumn])
140                         dj = fabs(dj);
141                    else
142                         dj = 1.0e50;
143                    weight[iColumn] = dj;
144                    sort[iColumn] = iColumn;
145               }
146               // sort
147               CoinSort_2(weight, weight + numberColumns, sort);
148               // and fix others
149               for (int iColumn = smallNumberColumns; iColumn < numberColumns; iColumn++) {
150                    int kColumn = sort[iColumn];
151                    double value = solution[kColumn];
152                    columnLower[kColumn] = value;
153                    columnUpper[kColumn] = value;
154               }
155          }
156     }
157     delete [] weight;
158     delete [] sort;
159     delete [] saveLower;
160     delete [] saveUpper;
161     model.primal(1);
162     return 0;
163}
Note: See TracBrowser for help on using the repository browser.