数値積分 - それを並列化する方法は？

Question

私はOpenCLを使い始めています。私はベクトルの追加の例を見て理解することができました。しかし、私は台形法について考えていました。これは[a、b]のx^2の積分計算のコード（C）です。

double f(double x) 
{ 
    return x*x; 
} 

double Simple_Trap(double a, double b) 
{ 
    double fA, fB; 
    fA = f(a); 
    fB = f(b); 
    return ((fA + fB) * (b-a))/2; 
} 

double Comp_Trap(double a, double b) 
{ 
    double Suma = 0; 
    double i = 0; 
    i = a + INC; 
    Suma += Simple_Trap(a,i); 
    while(i < b) 
    { 
     i+=INC; 
     Suma += Simple_Trap(i,i + INC); 
    } 
    return Suma; 
} 

質問は、台形法を使用した積分計算のためのカーネルを得る方法ですか？

---

だから、私はアイデアを考えていた：パーシャル[i]は=は（A +オフセット）を統合し、その後Patrick87を述べたようにパーシャルの合計を計算するためにカーネルを作ります。

しかし、これが最善の方法ですか？

Answer 1

これは私が思いついたものです。私はこのカーネルのエンドツーエンドテストを行うことはできませんでした。私はもう少し時間があるときに更新を行います。

comp_trapは、上記のコードに基づいて基本分割方法&を征服する方法です。 comp_trap_multiは、各作業項目がサブセクションを分割するようにして精度を向上させます。

各作業グループが結果を返すために倍精度の配列を割り当てる必要があります。これは、回避したいベクトル割り当てを減らすのに役立ちます。

問題がある場合は教えてください。

更新：ダブルのOpenCLで

2オプションであるため

1）は、ハード、この値は私のシステムに最適である64にワークグループのサイズをコード化し、必要があります）、float型にすべての二重の参照を変更しました実験的に決定される。ホストプログラムが最終的にはターゲットシステム上の最適値だけを使用するため、この値をローカルフロート配列で渡すよりもハードコーディングするほうが好きです。

3）はい、私はコンセプトをキャッチが、どのようにそれを実装するには？間違った計算（a1は間違っていた、今良いはずです）

/* 
numerical-integration.cl 
*/ 

float f(float x) 
{ 
    return x*x; 
} 

float simple_trap(float a, float b) 
{ 
    float fA, fB; 
    fA = f(a); 
    fB = f(b); 
    return ((fA + fB) * (b-a))/2; 
} 

__kernel void comp_trap(
    float a, 
    float b, 
    __global float* sums) 
{ 
/* 
- assumes 1D global and local work dimensions 
- each work unit will calculate 1/get_global_size of the total sum 
- the 0th work unit of each group then accumulates the sum for the 
group and stores it in __global * sums 
- memory allocation: sizeof(sums) = get_num_groups(0) * sizeof(float) 
- assumes local scratchpad size is at lease 8 bytes per work unit in the group 
ie sizeof(wiSums) = get_local_size(0) * sizeof(float) 
*/ 
    __local float wiSums[64]; 
    int l_id = get_local_id(0); 

    //cumpute range for this work item is: a1, b1 
    float a1 = a+((b-a)/get_global_size(0))*get_global_id(0); 
    float b1 = a1+(b-a)/get_global_size(0); 

    wiSums[l_id] = simple_trap(a1,b1); 

    barrier(CLK_LOCAL_MEM_FENCE); 

    int i; 
    if(l_id == 0){ 
     for(i=1;i<get_local_size(0);i++){ 
      wiSums[0] += wiSums[i]; 
     } 
     sums[get_group_id(0)] = wiSums[0]; 
    } 
} 

__kernel void comp_trap_multi(
    float a, 
    float b, 
    __global float* sums, 
    int divisions) 
{ 
/* 
- same as above, but each work unit further divides its range into 
'divisions' equal parts, yielding a more accurate result 
- work units still store only one sum in the local array, which is 
used later for the final group accumulation 
*/ 
    __local float wiSums[64]; 
    int l_id = get_local_id(0); 

    float a1 = a+((b-a)/get_global_size(0))*get_global_id(0); 
    float b1 = a1+(b-a)/get_global_size(0); 
    float range; 
    if(divisions > 0){ 
     range = (b1-a1)/divisions; 
    }else{ 
     range = (b1-a1); 
    } 

    int i; 
    wiSums[l_id] = 0; 
    for(i=0;i<divisions;i++){ 
     wiSums[l_id] += simple_trap(a1+range*i,a1+range*(i+1)); 
    } 

    barrier(CLK_LOCAL_MEM_FENCE); 

    if(l_id == 0){ 
     for(i=1;i<get_local_size(0);i++){ 
      wiSums[0] += wiSums[i]; 
     } 
     sums[get_group_id(0)] = wiSums[0]; 
    } 
} 

Answer 2

台形法は、リーマン和の処理をわずかに改善したものです。これを並行して実行するには、間隔をスレッドと同じ数のサブインターバルに分割する必要があります。次に、各スレッドがそのサブインターバルで関数を統合するようにします。最後に、前のフェーズで計算されたすべての積分に対してグローバルな合計を減らします。各ステージで使用するスレッドの数を試すことができます。