Robert Crovella は、あなたの問題は反復回数に関して与えられた停止規則であるとすでに指摘しています。
最小限の高度な二分法では、停止規則はターゲットの精度にも関連する可能性があります。以下に、C++の数値レシピ本で利用可能なバージョンを適用した CUDA の二分法のバージョンを提供します。これにより、目標精度も設定できます。
eigenvaluesおそらく、 CUDA SDK サンプルで利用されている二分カーネルを適応させることで、計算的により洗練された二分法を取得できます。
メソッドの新しいバージョンは、より正確なようです。いくつかの結果の下:
No target accuracy
1.571289062500
4.453613281250
6.504882812500
10.546875000000
13.171386718750
Target accuracy
1.570796326795
4.712388980385
7.853981633975
10.995574287564
14.137166941154
Actual roots
1.570796326794897
4.712388980384690
7.853981633974483
10.995574287564276
14.137166941154069
上記の本で利用可能なものによって、より良い初期ブラケティングを再び達成することができます。
ここにコードがあります
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <math_constants.h>
#define BLOCKSIZE 512
/*******************/
/* iDivUp FUNCTION */
/*******************/
int iDivUp(int a, int b) { return ((a % b) != 0) ? (a / b + 1) : (a / b); }
/********************/
/* CUDA ERROR CHECK */
/********************/
#define gpuErrchk(ans) { gpuAssert((ans), __FILE__, __LINE__); }
inline void gpuAssert(cudaError_t code, char *file, int line, bool abort=true)
{
if (code != cudaSuccess)
{
fprintf(stderr,"GPUassert: %s %s %d\n", cudaGetErrorString(code), file, line);
if (abort) exit(code);
}
}
/************************************/
/* FUNCTION TO SEARCH THE ROOTS FOR */
/************************************/
__host__ __device__ double f(double x)
{
//return ((5*sin(2*x))-(52*cos(2*x)))+50;
return cos(x);
}
/***************************************/
/* BISECTION KERNEL - ORIGINAL VERSION */
/***************************************/
__global__ void bisection(double *a, double *b, double *c, int N, int loopcnt)
{
int tid = blockIdx.x*blockDim.x+threadIdx.x;
int loops = 0;
if (tid < N)
while (loops < loopcnt){
c[tid] = (a[tid] + b[tid])/2;
if ((f(c[tid]) * f(a[tid])) < 0) b[tid] = c[tid];
else a[tid] = c[tid];
loops++;
}
}
/************************************************/
/* BISECTION KERNEL - NUMERICAL RECIPES VERSION */
/************************************************/
// --- Using bisection, return the root of a function func known to lie between x1 and x2.
// The root will be refined until its accuracy is xacc.
__global__ void bisection_NR(const double *d_x1, const double *d_x2, double *d_roots, const double xacc, const int loopcnt, const int N) {
// --- loopcnt is the maximum allowed number of bisections.
int tid = blockIdx.x*blockDim.x+threadIdx.x;
if (tid < N) {
double dx,xmid,rtb;
double f1=f(d_x1[tid]);
double fmid=f(d_x2[tid]);
if (f1*fmid >= 0.0) d_roots[tid] = CUDART_NAN;
rtb = f1 < 0.0 ? (dx=d_x2[tid]-d_x1[tid],d_x1[tid]) : (dx=d_x1[tid]-d_x2[tid],d_x2[tid]); // --- Orient the search so that f>0
for (int j=0;j<loopcnt;j++) { // --- lies at x+dx.
fmid=f(xmid=rtb+(dx *= 0.5)); // --- Bisection loop.
if (fmid <= 0.0) rtb=xmid;
if (abs(dx) < xacc || fmid == 0.0) { d_roots[tid]=rtb; return; }
}
d_roots[tid] = CUDART_NAN;
}
}
/*******/
/* INT */
/*******/
int main()
{
int loops=100000; // --- Number of bisection iterations to run
double x1=0.0, x2=10.0; // --- Minimum and maximum values of the search interval
double Deltax = 1.0; // --- Sampling step of the search interval
int N = (int)x2/(int)Deltax; // --- Number of search intervales
// --- Host-side memory allocations
double *host_a = (double*)malloc(N*sizeof(double));
double *host_b = (double*)malloc(N*sizeof(double));
double *host_c = (double*)malloc(N*sizeof(double));
// --- Device-side memory allocations
double *dev_a; gpuErrchk(cudaMalloc(&dev_a, N*sizeof(double)));
double *dev_b; gpuErrchk(cudaMalloc(&dev_b, N*sizeof(double)));
double *dev_c; gpuErrchk(cudaMalloc(&dev_c, N*sizeof(double)));
// --- Initialize vectors on host
int i=0;
while(i < N) {
if((f(x1)*f(x1+Deltax))>0) x1 = x1 + Deltax;
else {
host_a[i] = x1;
host_b[i] = x1+Deltax;
x1 = x1 + Deltax;
i++;
}
}
// --- Copy host vectors to device
gpuErrchk(cudaMemcpy(dev_a, host_a, N*sizeof(double), cudaMemcpyHostToDevice));
gpuErrchk(cudaMemcpy(dev_b, host_b, N*sizeof(double), cudaMemcpyHostToDevice));
bisection<<<iDivUp(N, BLOCKSIZE), BLOCKSIZE>>>(dev_a, dev_b, dev_c, loops, N);
gpuErrchk(cudaMemcpy(host_c, dev_c, N*sizeof(double), cudaMemcpyDeviceToHost));
for(i=0; i<N; i++) printf("%3.12f\n",host_c[i]);
printf("\n");
bisection_NR<<<iDivUp(N, BLOCKSIZE), BLOCKSIZE>>>(dev_a, dev_b, dev_c, 2.5e-13, loops, N);
gpuErrchk(cudaMemcpy(host_c, dev_c, N*sizeof(double), cudaMemcpyDeviceToHost));
for(i=0; i<N; i++) printf("%3.12f\n",host_c[i]);
// --- Release device memory
gpuErrchk(cudaFree(dev_a));
gpuErrchk(cudaFree(dev_b));
gpuErrchk(cudaFree(dev_c));
// --- Release host memory
free(host_a);
free(host_b);
free(host_c);
return 0;
}