26 #ifdef SOPLEX_WITH_GMP 27 static void GMPv_init(mpz_t* vect,
int dim)
29 for(
int i = 0; i < dim; i++ )
33 static void GMPv_clear(mpz_t* vect,
int dim)
35 for(
int i = 0; i < dim; i++ )
41 static std::ostream&
operator<<(std::ostream& os,
const mpz_t* number)
43 os << mpz_get_str(0, 10, *number);
53 static int Reconstruct(
VectorRational& resvec, mpz_t* xnum, mpz_t denom,
int dim,
const Rational& denomBoundSquared,
const DIdxSet* indexSet = 0)
59 assert(mpz_sgn(denom) > 0);
60 assert(denomBoundSquared > 0);
76 mpz_set_q(Dbound, denomBoundSquared.getMpqRef());
78 mpz_set(Dbound, denom);
81 mpz_sqrt(Dbound, Dbound);
83 MSG_DEBUG( std::cout <<
"reconstructing " << dim <<
" dimensional vector with denominator bound " << mpz_get_str(0, 10, Dbound) <<
"\n" );
88 if( mpz_cmp_ui(Dbound,16777216) < 0 )
89 mpz_set_ui(Dbound,16777216);
103 for(
int c = 0; (indexSet == 0 && c < dim) || (indexSet != 0 && c < indexSet->size()); c++ )
105 int j = (indexSet == 0 ? c : indexSet->index(c));
110 MSG_DEBUG( std::cout <<
" --> component " << j <<
" = " << &xnum[j] <<
" / denom\n" );
113 if( mpz_sgn(xnum[j]) != 0 )
121 mpz_divexact(tn,tn,temp);
122 mpz_divexact(td,td,temp);
124 if(mpz_cmp(td,Dbound)<=0)
128 mpq_set_num(resvec[j].getMpqRef_w(),tn);
129 mpq_set_den(resvec[j].getMpqRef_w(),td);
137 mpz_fdiv_q(a0,tn,td);
138 mpz_fdiv_r(temp,tn,td);
141 mpz_fdiv_q(ai,tn,td);
142 mpz_fdiv_r(temp,tn,td);
148 mpz_addmul(p[2],a0,ai);
156 if( mpz_cmp(q[2],Dbound) > 0 )
163 while( !done && mpz_cmp_ui(td,0) )
168 mpz_fdiv_q(ai, tn, td);
169 mpz_fdiv_r(temp, tn, td);
181 mpz_addmul(p[2], p[1], ai);
183 mpz_addmul(q[2], q[1], ai);
185 if( mpz_cmp(q[2], Dbound) > 0 )
189 MSG_DEBUG( std::cout <<
" --> convergent denominator = " << &q[2] <<
"\n" );
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);
199 if( mpz_cmp(gcd, Dbound) > 0 )
201 MSG_DEBUG( std::cout <<
"terminating with gcd " << &gcd <<
" exceeding Dbound " << &Dbound <<
"\n" );
229 #ifdef SOPLEX_WITH_GMP 239 GMPv_init(xnum, dim);
240 mpz_init_set_ui(denom, 1);
245 for(
int i = 0; i < dim; i++ )
246 mpz_lcm(denom, denom, mpq_denref(input[i].getMpqRef()));
248 for(
int i = 0; i < dim; i++ )
250 mpz_mul(xnum[i], denom, mpq_numref(input[i].getMpqRef()));
251 mpz_divexact(xnum[i], xnum[i], mpq_denref(input[i].getMpqRef()));
256 for(
int i = 0; i < indexSet->
size(); i++ )
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()));
263 for(
int i = 0; i < indexSet->
size(); i++ )
265 int k = indexSet->
index(i);
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()));
273 MSG_DEBUG( std::cout <<
"LCM = " << mpz_get_str(0, 10, denom) <<
"\n" );
276 rval = Reconstruct(input, xnum, denom, dim, denomBoundSquared, indexSet);
279 GMPv_clear(xnum, dim);
297 if( solution.hasPrimal() )
310 if( solution.hasPrimalray() )
312 buffer.
reDim((solution._primalray).dim());
313 solution.getPrimalray(buffer);
315 solution._primalray = buffer;
318 if( solution.hasDual() )
323 solution.
_dual = buffer;
325 buffer.
reDim((solution._redcost).dim());
326 solution.getRedcost(buffer);
328 solution._redcost = buffer;
331 if( solution.hasDualfarkas() )
333 buffer.
reDim((solution._dualfarkas).dim());
334 solution.getDualfarkas(buffer);
336 solution._dualfarkas = buffer;
bool getDual(VectorBase< R > &vector) const
gets the dual solution vector; returns true on success
void reDim(int newdim, const bool setZero=true)
Resets DVectorBase's dimension to newdim.
Dense vector.Class VectorBase provides dense linear algebra vectors. It does not provide memory manag...
bool getSlacks(VectorBase< R > &vector) const
gets the vector of slack values; returns true on success
Wrapper for GMP type mpq_class.We wrap mpq_class so that we can replace it by a double type if GMP is...
bool getPrimal(VectorBase< R > &vector) const
gets the primal solution vector; returns true on success
void spx_alloc(T &p, int n=1)
Allocate memory.
Rational reconstruction of solution vector.
bool reconstructVector(VectorRational &input, const Rational &denomBoundSquared, const DIdxSet *indexSet)
bool reconstructSol(SolRational &solution)
std::ostream & operator<<(std::ostream &s, const VectorBase< R > &vec)
Output operator.
int size() const
returns the number of used indices.
Dynamic index set.Class DIdxSet provides dynamic IdxSet in the sense, that no restrictions are posed ...
Debugging, floating point type and parameter definitions.
int dim() const
Dimension of vector.
Everything should be within this namespace.
VectorBase< Rational > VectorRational
int index(int n) const
access n 'th index.
void spx_free(T &p)
Release memory.