mirror/eigen
2
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-01-30 17:40:05 +08:00
Code Issues Packages Projects Releases Wiki Activity

Finishing touches on igamma/igammac for GPU. Tests now pass.

Browse Source
This commit is contained in:
Eugene Brevdo 2016-03-07 15:35:09 -08:00
parent 5707004d6b
commit 0bb5de05a1
2 changed files with 16 additions and 18 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

22
Eigen/src/Core/SpecialFunctions.h
View File

@ -520,14 +520,16 @@ struct igammac_impl {
Copyright 1985, 1987, 1992 by Stephen L. Moshier Copyright 1985, 1987, 1992 by Stephen L. Moshier
Direct inquiries to 30 Frost Street, Cambridge, MA 02140 Direct inquiries to 30 Frost Street, Cambridge, MA 02140
*/ */
using std::log;
const Scalar zero = 0; const Scalar zero = 0;
const Scalar one = 1; const Scalar one = 1;
const Scalar two = 2; const Scalar two = 2;
const Scalar machep = igamma_helper<Scalar>::machep(); const Scalar machep = igamma_helper<Scalar>::machep();
const Scalar maxlog = ::log(NumTraits<Scalar>::highest()); const Scalar maxlog = log(NumTraits<Scalar>::highest());
const Scalar big = igamma_helper<Scalar>::big(); const Scalar big = igamma_helper<Scalar>::big();
const Scalar biginv = 1 / big; const Scalar biginv = 1 / big;
const Scalar nan = NumTraits<Scalar>::quiet_NaN(); const Scalar nan = NumTraits<Scalar>::quiet_NaN();
const Scalar inf = NumTraits<Scalar>::infinity();
Scalar ans, ax, c, yc, r, t, y, z; Scalar ans, ax, c, yc, r, t, y, z;
Scalar pk, pkm1, pkm2, qk, qkm1, qkm2; Scalar pk, pkm1, pkm2, qk, qkm1, qkm2;
@ -541,11 +543,9 @@ struct igammac_impl {
return (one - igamma_impl<Scalar>::run(a, x)); return (one - igamma_impl<Scalar>::run(a, x));
} }
using std::isinf; if (x == inf) return zero; // std::isinf crashes on CUDA
if ((isinf)(x)) return zero;
/* Compute x**a * exp(-x) / gamma(a) */ /* Compute x**a * exp(-x) / gamma(a) */
using std::log;
ax = a * log(x) - x - lgamma_impl<Scalar>::run(a); ax = a * log(x) - x - lgamma_impl<Scalar>::run(a);
if (ax < -maxlog) { // underflow if (ax < -maxlog) { // underflow
return zero; return zero;
@ -564,7 +564,7 @@ struct igammac_impl {
ans = pkm1 / qkm1; ans = pkm1 / qkm1;
using std::abs; using std::abs;
do { while (true) {
c += one; c += one;
y += one; y += one;
z += two; z += two;
@ -588,7 +588,8 @@ struct igammac_impl {
qkm2 *= biginv; qkm2 *= biginv;
qkm1 *= biginv; qkm1 *= biginv;
} }
} while (t > machep); if (t <= machep) break;
}
return (ans * ax); return (ans * ax);
} }
@ -683,10 +684,11 @@ struct igamma_impl {
* k=0 | (a+k+1) * k=0 | (a+k+1)
* *
*/ */
using std::log;
const Scalar zero = 0; const Scalar zero = 0;
const Scalar one = 1; const Scalar one = 1;
const Scalar machep = igamma_helper<Scalar>::machep(); const Scalar machep = igamma_helper<Scalar>::machep();
const Scalar maxlog = ::log(NumTraits<Scalar>::highest()); const Scalar maxlog = log(NumTraits<Scalar>::highest());
const Scalar nan = NumTraits<Scalar>::quiet_NaN(); const Scalar nan = NumTraits<Scalar>::quiet_NaN();
double ans, ax, c, r; double ans, ax, c, r;
@ -702,7 +704,6 @@ struct igamma_impl {
} }
/* Compute x**a * exp(-x) / gamma(a) */ /* Compute x**a * exp(-x) / gamma(a) */
using std::log;
ax = a * log(x) - x - lgamma_impl<Scalar>::run(a); ax = a * log(x) - x - lgamma_impl<Scalar>::run(a);
if (ax < -maxlog) { if (ax < -maxlog) {
// underflow // underflow
@ -716,11 +717,12 @@ struct igamma_impl {
c = one; c = one;
ans = one; ans = one;
do { while (true) {
r += one; r += one;
c *= x/r; c *= x/r;
ans += c; ans += c;
} while (c/ans > machep); if (c/ans <= machep) break;
}
return (ans * ax / a); return (ans * ax / a);
} }

12
unsupported/test/cxx11_tensor_cuda.cu
View File

@ -632,7 +632,6 @@ void test_cuda_igamma()
Tensor<Scalar, 2> a(6, 6); Tensor<Scalar, 2> a(6, 6);
Tensor<Scalar, 2> x(6, 6); Tensor<Scalar, 2> x(6, 6);
Tensor<Scalar, 2> out(6, 6); Tensor<Scalar, 2> out(6, 6);
Tensor<Scalar, 2> expected_out(6, 6);
out.setZero(); out.setZero();
Scalar a_s[] = {Scalar(0), Scalar(1), Scalar(1.5), Scalar(4), Scalar(0.0001), Scalar(1000.5)}; Scalar a_s[] = {Scalar(0), Scalar(1), Scalar(1.5), Scalar(4), Scalar(0.0001), Scalar(1000.5)};
@ -641,7 +640,7 @@ void test_cuda_igamma()
for (int i = 0; i < 6; ++i) { for (int i = 0; i < 6; ++i) {
for (int j = 0; j < 6; ++j) { for (int j = 0; j < 6; ++j) {
a(i, j) = a_s[i]; a(i, j) = a_s[i];
x(i, j) = x_s[i]; x(i, j) = x_s[j];
} }
} }
@ -686,8 +685,7 @@ void test_cuda_igamma()
for (int i = 0; i < 6; ++i) { for (int i = 0; i < 6; ++i) {
for (int j = 0; j < 6; ++j) { for (int j = 0; j < 6; ++j) {
if ((std::isnan)(igamma_s[i][j])) { if ((std::isnan)(igamma_s[i][j])) {
printf("got: %g\n", out(i, j)); VERIFY((std::isnan)(out(i, j)));
//VERIFY((std::isnan)(out(i, j)));
} else { } else {
VERIFY_IS_APPROX(out(i, j), igamma_s[i][j]); VERIFY_IS_APPROX(out(i, j), igamma_s[i][j]);
} }
@ -701,7 +699,6 @@ void test_cuda_igammac()
Tensor<Scalar, 2> a(6, 6); Tensor<Scalar, 2> a(6, 6);
Tensor<Scalar, 2> x(6, 6); Tensor<Scalar, 2> x(6, 6);
Tensor<Scalar, 2> out(6, 6); Tensor<Scalar, 2> out(6, 6);
Tensor<Scalar, 2> expected_out(6, 6);
out.setZero(); out.setZero();
Scalar a_s[] = {Scalar(0), Scalar(1), Scalar(1.5), Scalar(4), Scalar(0.0001), Scalar(1000.5)}; Scalar a_s[] = {Scalar(0), Scalar(1), Scalar(1.5), Scalar(4), Scalar(0.0001), Scalar(1000.5)};
@ -710,7 +707,7 @@ void test_cuda_igammac()
for (int i = 0; i < 6; ++i) { for (int i = 0; i < 6; ++i) {
for (int j = 0; j < 6; ++j) { for (int j = 0; j < 6; ++j) {
a(i, j) = a_s[i]; a(i, j) = a_s[i];
x(i, j) = x_s[i]; x(i, j) = x_s[j];
} }
} }
@ -754,8 +751,7 @@ void test_cuda_igammac()
for (int i = 0; i < 6; ++i) { for (int i = 0; i < 6; ++i) {
for (int j = 0; j < 6; ++j) { for (int j = 0; j < 6; ++j) {
if ((std::isnan)(igammac_s[i][j])) { if ((std::isnan)(igammac_s[i][j])) {
printf("got: %g\n", out(i, j)); VERIFY((std::isnan)(out(i, j)));
//VERIFY((std::isnan)(out(i, j)));
} else { } else {
VERIFY_IS_APPROX(out(i, j), igammac_s[i][j]); VERIFY_IS_APPROX(out(i, j), igammac_s[i][j]);
} }