|
Go to the documentation of this file.
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);
138 mpz_fdiv_q(a0,tn,td);
139 mpz_fdiv_r(temp,tn,td);
142 mpz_fdiv_q(ai,tn,td);
143 mpz_fdiv_r(temp,tn,td);
149 mpz_addmul(p[2],a0,ai);
157 if( mpz_cmp(q[2],Dbound) > 0 )
164 while( !done && mpz_cmp_ui(td,0) )
169 mpz_fdiv_q(ai, tn, td);
170 mpz_fdiv_r(temp, tn, td);
182 mpz_addmul(p[2], p[1], ai);
184 mpz_addmul(q[2], q[1], ai);
186 if( mpz_cmp(q[2], Dbound) > 0 )
190 MSG_DEBUG( std::cout << " --> convergent denominator = " << &q[2] << "\n" );
194 mpq_set_num(resvec[j].getMpqRef_w(), p[1]);
195 mpq_set_den(resvec[j].getMpqRef_w(), q[1]);
196 mpq_canonicalize(resvec[j].getMpqRef_w());
197 mpz_gcd(temp, gcd, mpq_denref(resvec[j].getMpqRef()));
198 mpz_mul(gcd, gcd, temp);
200 if( mpz_cmp(gcd, Dbound) > 0 )
202 MSG_DEBUG( std::cout << "terminating with gcd " << &gcd << " exceeding Dbound " << &Dbound << "\n" );
230 #ifdef SOPLEX_WITH_GMP
240 GMPv_init(xnum, dim);
241 mpz_init_set_ui(denom, 1);
246 for( int i = 0; i < dim; i++ )
247 mpz_lcm(denom, denom, mpq_denref(input[i].getMpqRef()));
249 for( int i = 0; i < dim; i++ )
251 mpz_mul(xnum[i], denom, mpq_numref(input[i].getMpqRef()));
252 mpz_divexact(xnum[i], xnum[i], mpq_denref(input[i].getMpqRef()));
257 for( int i = 0; i < indexSet-> size(); i++ )
259 assert(indexSet-> index(i) >= 0);
260 assert(indexSet-> index(i) < input. dim());
261 mpz_lcm(denom, denom, mpq_denref(input[indexSet-> index(i)].getMpqRef()));
264 for( int i = 0; i < indexSet-> size(); i++ )
266 int k = indexSet-> index(i);
268 assert(k < input. dim());
269 mpz_mul(xnum[k], denom, mpq_numref(input[k].getMpqRef()));
270 mpz_divexact(xnum[k], xnum[k], mpq_denref(input[k].getMpqRef()));
274 MSG_DEBUG( std::cout << "LCM = " << mpz_get_str(0, 10, denom) << "\n" );
277 rval = Reconstruct(input, xnum, denom, dim, denomBoundSquared, indexSet);
280 GMPv_clear(xnum, dim);
311 if( solution.hasPrimalray() )
313 buffer. reDim((solution._primalray).dim());
314 solution.getPrimalray(buffer);
316 solution._primalray = buffer;
324 solution. _dual = buffer;
326 buffer. reDim((solution._redcost).dim());
327 solution.getRedcost(buffer);
329 solution._redcost = buffer;
332 if( solution.hasDualfarkas() )
334 buffer. reDim((solution._dualfarkas).dim());
335 solution.getDualfarkas(buffer);
337 solution._dualfarkas = buffer;
|