source: stable/1.13/Clp/examples/useVolume.cpp @ 1898

Last change on this file since 1898 was 1561, checked in by stefan, 10 years ago

make check for sample dir work without including ClpConfig?.h

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