Scippy

SoPlex

Sequential object-oriented simPlex

ratrecon.cpp
Go to the documentation of this file.
1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
2 /* */
3 /* This file is part of the class library */
4 /* SoPlex --- the Sequential object-oriented simPlex. */
5 /* */
6 /* Copyright (C) 1996-2017 Konrad-Zuse-Zentrum */
7 /* fuer Informationstechnik Berlin */
8 /* */
9 /* SoPlex is distributed under the terms of the ZIB Academic Licence. */
10 /* */
11 /* You should have received a copy of the ZIB Academic License */
12 /* along with SoPlex; see the file COPYING. If not email to soplex@zib.de. */
13 /* */
14 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
15 
16 #ifndef SOPLEX_LEGACY
17 
18 #include <iostream>
19 #include <assert.h>
20 
21 #include "spxdefines.h"
22 #include "ratrecon.h"
23 
24 namespace soplex
25 {
26 #ifdef SOPLEX_WITH_GMP
27  static void GMPv_init(mpz_t* vect, int dim)
28  {
29  for( int i = 0; i < dim; i++ )
30  mpz_init(vect[i]);
31  }
32 
33  static void GMPv_clear(mpz_t* vect, int dim)
34  {
35  for( int i = 0; i < dim; i++ )
36  mpz_clear(vect[i]);
37  }
38 
39 #ifdef SOPLEX_DEBUG
40  /** print integer to stream */
41  static std::ostream& operator<<(std::ostream& os, const mpz_t* number)
42  {
43  os << mpz_get_str(0, 10, *number);
44  return os;
45  }
46 #endif
47 
48  /** this reconstruction routine will set x equal to the mpq vector where each component is the best rational
49  * approximation of xnum / denom with where the GCD of denominators of x is at most Dbound; it will return true on
50  * success and false if more accuracy is required: specifically if componentwise rational reconstruction does not
51  * produce such a vector
52  */
53  static int Reconstruct(VectorRational& resvec, mpz_t* xnum, mpz_t denom, int dim, const Rational& denomBoundSquared, const DIdxSet* indexSet = 0)
54  {
55  bool rval = true;
56  int done = 0;
57 
58  /* denominator must be positive */
59  assert(mpz_sgn(denom) > 0);
60  assert(denomBoundSquared > 0);
61 
62  mpz_t temp;
63  mpz_t td;
64  mpz_t tn;
65  mpz_t Dbound;
66  mpz_t gcd;
67 
68  mpz_init(gcd); /* stores the gcd of denominators; abort if too big */
69  mpz_set_ui(gcd, 1);
70  mpz_init(temp);
71  mpz_init(td);
72  mpz_init(tn);
73  mpz_init(Dbound);
74 
75 #if 1
76  mpz_set_q(Dbound, denomBoundSquared.getMpqRef()); /* this is the working bound on the denominator size */
77 #else
78  mpz_set(Dbound, denom); /* this is the working bound on the denominator size */
79 #endif
80 
81  mpz_sqrt(Dbound, Dbound);
82 
83  MSG_DEBUG( std::cout << "reconstructing " << dim << " dimensional vector with denominator bound " << mpz_get_str(0, 10, Dbound) << "\n" );
84 
85  /* if Dbound is below 2^24 increase it to this value, this avoids changing input vectors that have low denominator
86  * because they are floating point representable
87  */
88  if( mpz_cmp_ui(Dbound,16777216) < 0 )
89  mpz_set_ui(Dbound,16777216);
90 
91  /* The following represent a_i, the cont frac representation and p_i/q_i, the convergents */
92  mpz_t a0;
93  mpz_t ai;
94  mpz_init(a0);
95  mpz_init(ai);
96 
97  /* here we use p[2]=pk, p[1]=pk-1,p[0]=pk-2 and same for q */
98  mpz_t p[3];
99  GMPv_init(p, 3);
100  mpz_t q[3];
101  GMPv_init(q, 3);
102 
103  for( int c = 0; (indexSet == 0 && c < dim) || (indexSet != 0 && c < indexSet->size()); c++ )
104  {
105  int j = (indexSet == 0 ? c : indexSet->index(c));
106 
107  assert(j >= 0);
108  assert(j < dim);
109 
110  MSG_DEBUG( std::cout << " --> component " << j << " = " << &xnum[j] << " / denom\n" );
111 
112  /* if xnum =0 , then just leave x[j] as zero */
113  if( mpz_sgn(xnum[j]) != 0 )
114  {
115  /* setup n and d for computing a_i the cont. frac. rep */
116  mpz_set(tn,xnum[j]);
117  mpz_set(td,denom);
118 
119  /* divide tn and td by gcd */
120  mpz_gcd(temp,tn,td);
121  mpz_divexact(tn,tn,temp);
122  mpz_divexact(td,td,temp);
123 
124  if(mpz_cmp(td,Dbound)<=0)
125  {
126  MSG_DEBUG( std::cout << "marker 1\n" );
127 
128  mpq_set_num(resvec[j].getMpqRef_w(),tn);
129  mpq_set_den(resvec[j].getMpqRef_w(),td);
130  }
131  else
132  {
133  MSG_DEBUG( std::cout << "marker 2\n" );
134 
135  mpz_set_ui(temp,1);
136 
137  mpz_fdiv_q(a0,tn,td);
138  mpz_fdiv_r(temp,tn,td);
139  mpz_set(tn,td);
140  mpz_set(td,temp);
141  mpz_fdiv_q(ai,tn,td);
142  mpz_fdiv_r(temp,tn,td);
143  mpz_set(tn,td);
144  mpz_set(td,temp);
145 
146  mpz_set(p[1],a0);
147  mpz_set_ui(p[2],1);
148  mpz_addmul(p[2],a0,ai);
149 
150  mpz_set_ui(q[1],1);
151  mpz_set(q[2],ai);
152 
153  done = 0;
154 
155  /* if q is already big, skip loop */
156  if( mpz_cmp(q[2],Dbound) > 0 )
157  {
158  MSG_DEBUG( std::cout << "marker 3\n" );
159  done = 1;
160  }
161 
162  int cfcnt = 2;
163  while( !done && mpz_cmp_ui(td,0) )
164  {
165  /* update everything: compute next ai, then update convergents */
166 
167  /* update ai */
168  mpz_fdiv_q(ai, tn, td);
169  mpz_fdiv_r(temp, tn, td);
170  mpz_set(tn, td);
171  mpz_set(td, temp);
172 
173  /* shift p,q */
174  mpz_set(q[0], q[1]);
175  mpz_set(q[1], q[2]);
176  mpz_set(p[0], p[1]);
177  mpz_set(p[1], p[2]);
178 
179  /* compute next p,q */
180  mpz_set(p[2], p[0]);
181  mpz_addmul(p[2], p[1], ai);
182  mpz_set(q[2], q[0]);
183  mpz_addmul(q[2], q[1], ai);
184 
185  if( mpz_cmp(q[2], Dbound) > 0 )
186  done = 1;
187  cfcnt++;
188 
189  MSG_DEBUG( std::cout << " --> convergent denominator = " << &q[2] << "\n" );
190  }
191 
192  /* Assign the values */
193  mpq_set_num(resvec[j].getMpqRef_w(), p[1]);
194  mpq_set_den(resvec[j].getMpqRef_w(), q[1]);
195  mpq_canonicalize(resvec[j].getMpqRef_w());
196  mpz_gcd(temp, gcd, mpq_denref(resvec[j].getMpqRef()));
197  mpz_mul(gcd, gcd, temp);
198 
199  if( mpz_cmp(gcd, Dbound) > 0 )
200  {
201  MSG_DEBUG( std::cout << "terminating with gcd " << &gcd << " exceeding Dbound " << &Dbound << "\n" );
202  rval = false;
203  goto CLEANUP;
204  }
205  }
206  }
207  }
208 
209  CLEANUP:
210  GMPv_clear(q, 3);
211  GMPv_clear(p, 3);
212  mpz_clear(td);
213  mpz_clear(tn);
214  mpz_clear(a0);
215  mpz_clear(ai);
216  mpz_clear(temp);
217  mpz_clear(Dbound);
218  mpz_clear(gcd);
219 
220  return rval;
221  }
222 #endif
223 
224 
225 
226  /** reconstruct a rational vector */
227  bool reconstructVector(VectorRational& input, const Rational& denomBoundSquared, const DIdxSet* indexSet)
228  {
229 #ifdef SOPLEX_WITH_GMP
230  mpz_t* xnum = 0; /* numerator of input vector */
231  mpz_t denom; /* common denominator of input vector */
232  int rval = true;
233  int dim;
234 
235  dim = input.dim();
236 
237  /* convert vector to mpz format */
238  spx_alloc(xnum, dim);
239  GMPv_init(xnum, dim);
240  mpz_init_set_ui(denom, 1);
241 
242  /* find common denominator */
243  if( indexSet == 0 )
244  {
245  for( int i = 0; i < dim; i++ )
246  mpz_lcm(denom, denom, mpq_denref(input[i].getMpqRef()));
247 
248  for( int i = 0; i < dim; i++ )
249  {
250  mpz_mul(xnum[i], denom, mpq_numref(input[i].getMpqRef()));
251  mpz_divexact(xnum[i], xnum[i], mpq_denref(input[i].getMpqRef()));
252  }
253  }
254  else
255  {
256  for( int i = 0; i < indexSet->size(); i++ )
257  {
258  assert(indexSet->index(i) >= 0);
259  assert(indexSet->index(i) < input.dim());
260  mpz_lcm(denom, denom, mpq_denref(input[indexSet->index(i)].getMpqRef()));
261  }
262 
263  for( int i = 0; i < indexSet->size(); i++ )
264  {
265  int k = indexSet->index(i);
266  assert(k >= 0);
267  assert(k < input.dim());
268  mpz_mul(xnum[k], denom, mpq_numref(input[k].getMpqRef()));
269  mpz_divexact(xnum[k], xnum[k], mpq_denref(input[k].getMpqRef()));
270  }
271  }
272 
273  MSG_DEBUG( std::cout << "LCM = " << mpz_get_str(0, 10, denom) << "\n" );
274 
275  /* reconstruct */
276  rval = Reconstruct(input, xnum, denom, dim, denomBoundSquared, indexSet);
277 
278  mpz_clear(denom);
279  GMPv_clear(xnum, dim);
280  spx_free(xnum);
281 
282  return rval;
283 #else
284  return false;
285 #endif
286  }
287 
288 
289 
290  /** reconstruct a rational solution */
291  /**@todo make this a method of class SoPlex */
292  bool reconstructSol(SolRational& solution)
293  {
294 #if 0
295  DVectorRational buffer;
296 
297  if( solution.hasPrimal() )
298  {
299  buffer.reDim((solution._primal).dim());
300  solution.getPrimal(buffer);
301  reconstructVector(buffer);
302  solution._primal = buffer;
303 
304  buffer.reDim((solution._slacks).dim());
305  solution.getSlacks(buffer);
306  reconstructVector(buffer);
307  solution._slacks = buffer;
308  }
309 
310  if( solution.hasPrimalray() )
311  {
312  buffer.reDim((solution._primalray).dim());
313  solution.getPrimalray(buffer);
314  reconstructVector(buffer);
315  solution._primalray = buffer;
316  }
317 
318  if( solution.hasDual() )
319  {
320  buffer.reDim((solution._dual).dim());
321  solution.getDual(buffer);
322  reconstructVector(buffer);
323  solution._dual = buffer;
324 
325  buffer.reDim((solution._redcost).dim());
326  solution.getRedcost(buffer);
327  reconstructVector(buffer);
328  solution._redcost = buffer;
329  }
330 
331  if( solution.hasDualfarkas() )
332  {
333  buffer.reDim((solution._dualfarkas).dim());
334  solution.getDualfarkas(buffer);
335  reconstructVector(buffer);
336  solution._dualfarkas = buffer;
337  }
338 #endif
339  return true;
340  }
341 } // namespace soplex
342 
343 #endif
DVectorBase< R > _slacks
Definition: solbase.h:219
bool getDual(VectorBase< R > &vector) const
gets the dual solution vector; returns true on success
Definition: solbase.h:94
void reDim(int newdim, const bool setZero=true)
Resets DVectorBase&#39;s dimension to newdim.
Definition: dvectorbase.h:249
Dense vector.Class VectorBase provides dense linear algebra vectors. It does not provide memory manag...
Definition: dsvectorbase.h:28
bool getSlacks(VectorBase< R > &vector) const
gets the vector of slack values; returns true on success
Definition: solbase.h:65
Wrapper for GMP type mpq_class.We wrap mpq_class so that we can replace it by a double type if GMP is...
Definition: rational.h:45
bool getPrimal(VectorBase< R > &vector) const
gets the primal solution vector; returns true on success
Definition: solbase.h:57
void spx_alloc(T &p, int n=1)
Allocate memory.
Definition: spxalloc.h:48
Rational reconstruction of solution vector.
bool reconstructVector(VectorRational &input, const Rational &denomBoundSquared, const DIdxSet *indexSet)
Definition: ratrecon.cpp:227
DVectorBase< R > _primal
Definition: solbase.h:218
bool reconstructSol(SolRational &solution)
Definition: ratrecon.cpp:292
#define MSG_DEBUG(x)
Definition: spxdefines.h:132
std::ostream & operator<<(std::ostream &s, const VectorBase< R > &vec)
Output operator.
Definition: basevectors.h:1159
int size() const
returns the number of used indices.
Definition: idxset.h:124
Dynamic index set.Class DIdxSet provides dynamic IdxSet in the sense, that no restrictions are posed ...
Definition: didxset.h:42
Debugging, floating point type and parameter definitions.
int dim() const
Dimension of vector.
Definition: vectorbase.h:215
Everything should be within this namespace.
VectorBase< Rational > VectorRational
Definition: vector.h:31
int index(int n) const
access n &#39;th index.
Definition: idxset.h:118
DVectorBase< R > _dual
Definition: solbase.h:221
void spx_free(T &p)
Release memory.
Definition: spxalloc.h:109