summaryrefslogtreecommitdiff
path: root/Build/source/utils/asymptote/cudareflect/simpson.cuh
diff options
context:
space:
mode:
Diffstat (limited to 'Build/source/utils/asymptote/cudareflect/simpson.cuh')
-rw-r--r--Build/source/utils/asymptote/cudareflect/simpson.cuh103
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;
+ }
+ }
+}
+