diff options
Diffstat (limited to 'Build/source/utils/asymptote/cudareflect/simpson.cuh')
-rw-r--r-- | Build/source/utils/asymptote/cudareflect/simpson.cuh | 103 |
1 files changed, 103 insertions, 0 deletions
diff --git a/Build/source/utils/asymptote/cudareflect/simpson.cuh b/Build/source/utils/asymptote/cudareflect/simpson.cuh new file mode 100644 index 00000000000..99fe8b7aad6 --- /dev/null +++ b/Build/source/utils/asymptote/cudareflect/simpson.cuh @@ -0,0 +1,103 @@ +/* +* John C. Bowman and Supakorn "Jamie" Rassameemasmuang <jamievlin@outlook.com> +* University of Alberta +* CUDA Adaptive Simpson integration +*/ + +#include <cuda.h> + +#include "utils.cuh" + +// Compute a numerical approximation to an integral via adaptive Simpson's Rule +// This routine ignores underflow. + +__device__ constexpr float sixth=1.0/6.0; +__device__ constexpr int depth=10; +__device__ constexpr int acc = 1.0/256; + +template<typename T> +struct TABLE { + bool left; // left interval? + float dat; + T psum, f1t, f2t, f3t, estr; +}; + +template<typename T=glm::vec3, typename TUtil=Vec3Utility, typename Tf> +__device__ inline T +simpson(Tf f, // Function to be integrated. + float a, float b, // Lower, upper limits of integration. + float acc) // Desired relative accuracy of integral. + // Try to make |error| <= acc*abs(integral). +{ + T integral,diff,area,estl,estr,est,fv0,fv1,fv2,fv3,fv4; + float dx; + TABLE<T> table[depth],*p,*pstop; + + p=table; + pstop=table+depth-1; + p->left=true; + p->psum=TUtil::init(); + float alpha=a; + float da=b-a; + fv0=f(alpha); + fv2=f(alpha+0.5f*da); + fv4=f(alpha+da); + float wt=sixth*da; + est=wt*(fv0+4.0f*fv2+fv4); + area=est; + float acc2=acc*acc; + + // Have estimate est of integral on (alpha, alpha+da). + // Bisect and compute estimates on left and right half intervals. + // integral is the best value for the integral. + + for(;;) { + dx=0.5f*da; + float arg=alpha+0.5f*dx; + fv1=f(arg); + fv3=f(arg+dx); + wt=sixth*dx; + estl=wt*(fv0+4.0f*fv1+fv2); + estr=wt*(fv2+4.0f*fv3+fv4); + integral=estl+estr; + diff=est-integral; + area -= diff; + + if(p >= pstop || (TUtil::abs2(diff) <= acc2*TUtil::abs2(area))) { + // Accept approximate integral. + // If it was a right interval, add results to finish at this level. + // If it was a left interval, process right interval. + + for(;;) { + if(p->left == false) { // process right-half interval + alpha += da; + p->left=true; + p->psum=integral; + fv0=p->f1t; + fv2=p->f2t; + fv4=p->f3t; + da=p->dat; + est=p->estr; + break; + } + integral += p->psum; + if(--p <= table) return integral; + } + + } else { + // Raise level and store information for processing right-half interval. + ++p; + da=dx; + est=estl; + p->left=false; + p->f1t=fv2; + p->f2t=fv3; + p->f3t=fv4; + p->dat=dx; + p->estr=estr; + fv4=fv2; + fv2=fv1; + } + } +} + |