source: trunk/Clp/examples/useVolume.cpp

Last change on this file was 1935, checked in by stefan, 5 years ago

fix compiler warnings; fix path to sample instance in minimum

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 9.6 KB
Line 
1/* $Id: useVolume.cpp 1935 2013-04-08 19:16:07Z 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 "ClpFactorization.hpp"
8#include "VolVolume.hpp"
9
10//#############################################################################
11
12class lpHook : public VOL_user_hooks {
13private:
14     lpHook(const lpHook&);
15     lpHook& operator= (const lpHook&);
16private:
17     /// Pointer to dense vector of structural variable upper bounds
18     double  *colupper_;
19     /// Pointer to dense vector of structural variable lower bounds
20     double  *collower_;
21     /// Pointer to dense vector of objective coefficients
22     double  *objcoeffs_;
23     /// Pointer to dense vector of right hand sides
24     double  *rhs_;
25     /// Pointer to dense vector of senses
26     char    *sense_;
27
28     /// The problem matrix in a row ordered form
29     CoinPackedMatrix rowMatrix_;
30     /// The problem matrix in a column ordered form
31     CoinPackedMatrix colMatrix_;
32
33public:
34     lpHook(const double* clb, const double* cub, const double* obj,
35            const double* rhs, const char* sense, const CoinPackedMatrix& mat);
36     virtual ~lpHook();
37
38public:
39     // for all hooks: return value of -1 means that volume should quit
40     /** compute reduced costs
41         @param u (IN) the dual variables
42         @param rc (OUT) the reduced cost with respect to the dual values
43     */
44     virtual int compute_rc(const VOL_dvector& u, VOL_dvector& rc);
45
46     /** Solve the subproblem for the subgradient step.
47         @param dual (IN) the dual variables
48         @param rc (IN) the reduced cost with respect to the dual values
49         @param lcost (OUT) the lagrangean cost with respect to the dual values
50         @param x (OUT) the primal result of solving the subproblem
51         @param v (OUT) b-Ax for the relaxed constraints
52         @param pcost (OUT) the primal objective value of <code>x</code>
53     */
54     virtual int solve_subproblem(const VOL_dvector& dual, const VOL_dvector& rc,
55                                  double& lcost, VOL_dvector& x, VOL_dvector& v,
56                                  double& pcost);
57     /** Starting from the primal vector x, run a heuristic to produce
58         an integer solution
59         @param x (IN) the primal vector
60         @param heur_val (OUT) the value of the integer solution (return
61         <code>DBL_MAX</code> here if no feas sol was found
62     */
63     virtual int heuristics(const VOL_problem& p,
64                            const VOL_dvector& x, double& heur_val) {
65          return 0;
66     }
67};
68
69//#############################################################################
70
71lpHook::lpHook(const double* clb, const double* cub, const double* obj,
72               const double* rhs, const char* sense,
73               const CoinPackedMatrix& mat)
74{
75     const int colnum = mat.getNumCols();
76     const int rownum = mat.getNumRows();
77
78     colupper_ = new double[colnum];
79     collower_ = new double[colnum];
80     objcoeffs_ = new double[colnum];
81     rhs_ = new double[rownum];
82     sense_ = new char[rownum];
83
84     std::copy(clb, clb + colnum, collower_);
85     std::copy(cub, cub + colnum, colupper_);
86     std::copy(obj, obj + colnum, objcoeffs_);
87     std::copy(rhs, rhs + rownum, rhs_);
88     std::copy(sense, sense + rownum, sense_);
89
90     if (mat.isColOrdered()) {
91          colMatrix_.copyOf(mat);
92          rowMatrix_.reverseOrderedCopyOf(mat);
93     } else {
94          rowMatrix_.copyOf(mat);
95          colMatrix_.reverseOrderedCopyOf(mat);
96     }
97}
98
99//-----------------------------------------------------------------------------
100
101lpHook::~lpHook()
102{
103     delete[] colupper_;
104     delete[] collower_;
105     delete[] objcoeffs_;
106     delete[] rhs_;
107     delete[] sense_;
108}
109
110//#############################################################################
111
112int
113lpHook::compute_rc(const VOL_dvector& u, VOL_dvector& rc)
114{
115     rowMatrix_.transposeTimes(u.v, rc.v);
116     const int psize = rowMatrix_.getNumCols();
117
118     for (int i = 0; i < psize; ++i)
119          rc[i] = objcoeffs_[i] - rc[i];
120     return 0;
121}
122
123//-----------------------------------------------------------------------------
124
125int
126lpHook::solve_subproblem(const VOL_dvector& dual, const VOL_dvector& rc,
127                         double& lcost, VOL_dvector& x, VOL_dvector& v,
128                         double& pcost)
129{
130     int i;
131     const int psize = x.size();
132     const int dsize = v.size();
133
134     // compute the lagrangean solution corresponding to the reduced costs
135     for (i = 0; i < psize; ++i)
136          x[i] = (rc[i] >= 0.0) ? collower_[i] : colupper_[i];
137
138     // compute the lagrangean value (rhs*dual + primal*rc)
139     lcost = 0;
140     for (i = 0; i < dsize; ++i)
141          lcost += rhs_[i] * dual[i];
142     for (i = 0; i < psize; ++i)
143          lcost += x[i] * rc[i];
144
145     // compute the rhs - lhs
146     colMatrix_.times(x.v, v.v);
147     for (i = 0; i < dsize; ++i)
148          v[i] = rhs_[i] - v[i];
149
150     // compute the lagrangean primal objective
151     pcost = 0;
152     for (i = 0; i < psize; ++i)
153          pcost += x[i] * objcoeffs_[i];
154
155     return 0;
156}
157
158//#############################################################################
159
160int main(int argc, const char *argv[])
161{
162     ClpSimplex  model;
163     int status;
164
165     if (argc < 2) {
166#if defined(SAMPLEDIR)
167          status = model.readMps(SAMPLEDIR "/p0033.mps", true);
168#else
169          fprintf(stderr, "Do not know where to find sample MPS files.\n");
170          exit(1);
171#endif
172     } else
173          status = model.readMps(argv[1], true);
174     if( status != 0 )
175     {
176        printf("Error %d reading MPS file\n", status);
177        return status;
178     }
179     /*
180       This driver uses volume algorithm
181       then does dual - after adjusting costs
182       then solves real problem
183     */
184
185     // do volume for a bit
186     VOL_problem volprob;
187     const CoinPackedMatrix* mat = model.matrix();
188     const int psize = mat->getNumCols();
189     const int dsize = mat->getNumRows();
190     char * sense = new char[dsize];
191     double * rhs = new double[dsize];
192     const double * rowLower = model.rowLower();
193     const double * rowUpper = model.rowUpper();
194     // Set the lb/ub on the duals
195     volprob.dsize = dsize;
196     volprob.psize = psize;
197     volprob.dual_lb.allocate(dsize);
198     volprob.dual_ub.allocate(dsize);
199     volprob.dsol.allocate(dsize);
200     int i;
201     for (i = 0; i < dsize; ++i) {
202          if (rowUpper[i] == rowLower[i]) {
203               // 'E':
204               volprob.dual_lb[i] = -1.0e31;
205               volprob.dual_ub[i] = 1.0e31;
206               rhs[i] = rowUpper[i];
207               sense[i] = 'E';
208          } else if (rowLower[i] < -0.99e10 && rowUpper[i] < 0.99e10) {
209               // 'L':
210               volprob.dual_lb[i] = -1.0e31;
211               volprob.dual_ub[i] = 0.0;
212               rhs[i] = rowUpper[i];
213               sense[i] = 'L';
214          } else if (rowLower[i] > -0.99e10 && rowUpper[i] > 0.99e10) {
215               // 'G':
216               volprob.dual_lb[i] = 0.0;
217               volprob.dual_ub[i] = 1.0e31;
218               rhs[i] = rowLower[i];
219               sense[i] = 'G';
220          } else {
221               printf("Volume Algorithm can't work if there is a non ELG row\n");
222               abort();
223          }
224     }
225     // Can't use read_param as private
226     // anyway I want automatic use - so maybe this is problem
227#if 0
228     FILE* infile = fopen("parameters", "r");
229     if (!infile) {
230          printf("Failure to open parameter file\n");
231     } else {
232          volprob.read_params("parameters");
233     }
234#endif
235#if 0
236     // should save and restore bounds
237     model.tightenPrimalBounds();
238#else
239     double * colUpper = model.columnUpper();
240     for (i = 0; i < psize; i++)
241          colUpper[i] = 1.0;
242#endif
243     lpHook myHook(model.getColLower(), model.getColUpper(),
244                   model.getObjCoefficients(),
245                   rhs, sense, *mat);
246     // move duals
247     double * pi = model.dualRowSolution();
248     memcpy(volprob.dsol.v, pi, dsize * sizeof(double));
249     volprob.solve(myHook,  false /* not warmstart */);
250     // For now stop as not doing any good
251     exit(77);
252     // create objectives
253     int numberRows = model.numberRows();
254     int numberColumns = model.numberColumns();
255     memcpy(pi, volprob.dsol.v, numberRows * sizeof(double));
256#define MODIFYCOSTS
257#ifdef MODIFYCOSTS
258     double * saveObj = new double[numberColumns];
259     memcpy(saveObj, model.objective(), numberColumns * sizeof(double));
260     memcpy(model.dualColumnSolution(), model.objective(),
261            numberColumns * sizeof(double));
262     model.clpMatrix()->transposeTimes(-1.0, pi, model.dualColumnSolution());
263     memcpy(model.objective(), model.dualColumnSolution(),
264            numberColumns * sizeof(double));
265     const double * rowsol = model.primalRowSolution();
266     //const double * rowLower = model.rowLower();
267     //const double * rowUpper = model.rowUpper();
268     double offset = 0.0;
269     for (i = 0; i < numberRows; i++) {
270          offset += pi[i] * rowsol[i];
271     }
272     double value2;
273     model.getDblParam(ClpObjOffset, value2);
274     printf("Offset %g %g\n", offset, value2);
275     model.setRowObjective(pi);
276     // zero out pi
277     memset(pi, 0, numberRows * sizeof(double));
278#endif
279     // Could put some in basis - only partially tested
280     model.allSlackBasis();
281     model.factorization()->maximumPivots(1000);
282     //model.setLogLevel(63);
283     // solve
284     model.dual(1);
285     //model.primal(1);
286#ifdef MODIFYCOSTS
287     memcpy(model.objective(), saveObj, numberColumns * sizeof(double));
288     // zero out pi
289     memset(pi, 0, numberRows * sizeof(double));
290     model.setRowObjective(pi);
291     delete [] saveObj;
292     model.primal();
293#endif
294
295     return 0;
296}
Note: See TracBrowser for help on using the repository browser.