BookGPU/Chapters/chapter11/code1.cu

   1
   2 template<typename Tx, typename Ty>
   3 __global__ void CalculateCoefficientsKnots( Tx *u, Ty *v, double *b, double *c,
   4   double *t, double *alpha, double *beta, double *gamma, int N )
   5 {
   6   int tid = threadIdx.x + blockIdx.x * blockDim.x;
   7   int s = tid*2;
   8   while(tid<=(N-2))
   9   {
  10     // decide whether an additional knot is necessary
  11     if(fabs(c[tid]+c[tid+1]- 2*b[tid])<=0.1e-5) // tolerance
  12     {
  13       //no additional knot
  14       h[s]=h[s+1]=u[tid];
  15       alpha[s]=alpha[s+1]=v[tid];
  16       beta[s]=beta[s+1]=c[tid];
  17       gamma[s]=gamma[s+1]=(c[tid+1]-c[tid])/(2*(fmax(1e-10,u[tid+1]-u[tid])));
  18     }
  19     else {
  20       //adding a knot
  21       h[s]=u[tid];
  22       //determine the position of the knot
  23       if((c[tid+1] - b[tid])*(c[tid] - b[tid])<0)
  24         h[s+1]=u[tid+1] + (c[tid] - b[tid])*(fmax(1e-10,u[tid+1]-u[tid]))/fmax(1e-10,(c[tid+1] - c[tid]));
  25       else
  26         h[s+1]=0.5*(u[tid+1] + u[tid]);
  27       //calculate coefficients
  28       double dtemp = (2*b[tid] - c[tid+1])+((c[tid+1] - c[tid])*(h[s+1] - u[tid]))/fmax(1e-10,(u[tid+1] - u[tid]));
  29       alpha[s]=v[tid];   beta[s]=c[tid];
  30       gamma[s]=(dtemp - c[tid])/(2*fmax(1e-10,(h[s+1] - u[tid])));
  31       alpha[s+1]=v[tid] + c[tid]*(h[s+1] - u[tid]) +
  32                 (dtemp - c[tid])*(h[s+1] - u[tid])/2;
  33       gamma[s+1]=(c[tid+1] - dtemp)/(2*fmax(1e-10,(u[tid+1] - h[s+1])));
  34       beta[s+1]=dtemp;
  35     }
  36     tid += blockDim.x * gridDim.x;   s = tid*2;
  37   }
  38   __syncthreads();
  39   // Select a single thread  to perform the last operation
  40   if((threadIdx.x  ) == 0)  {
  41     s = (N-1) * 2;   h[s]=u[N-1];
  42   }
  43   __syncthreads();
  44 }