3

元の質問

不均一なノード ポイントで補間を実行する次のカーネルがあり、最適化したいと考えています。

__global__ void interpolation(cufftDoubleComplex *Uj, double *points, cufftDoubleComplex *result, int N, int M)
{
    int i = threadIdx.x + blockDim.x * blockIdx.x;

    int PP;
    double P;
    const double alfa=(2.-1./cc)*pi_double-0.01;
    double phi_cap_s;
    cufftDoubleComplex temp;

    double cc_points=cc*points[i];
    double r_cc_points=rint(cc*points[i]);

    temp = make_cuDoubleComplex(0.,0.);

    if(i<M) {   
        for(int m=0; m<(2*K+1); m++) {
            P = (K*K-(cc_points-(r_cc_points+m-K))*(cc_points-(r_cc_points+m-K)));

            if(P>0.)  phi_cap_s = (1./pi_double)*((sinh(alfa*sqrt(P)))/sqrt(P));  
            if(P<0.)  phi_cap_s = (1./pi_double)*((sin(alfa*sqrt(-P)))/sqrt(-P));   
            if(P==0.) phi_cap_s = alfa/pi_double;        

            PP = modulo((r_cc_points + m -K ),(cc*N)); 
            temp.x = temp.x+phi_cap_s*Uj[PP].x; 
            temp.y = temp.y+phi_cap_s*Uj[PP].y; 
        } 

        result[i] = temp; 
    }
}

K と cc は定数で、points にはノードが含まれ、Uj には補間される値が含まれます。modulo は、基本的に % として機能する関数ですが、負の値に適切に拡張されます。特定の配置では、カーネル呼び出しに 2.3ms かかります。最も高価な部品は

            if(P>0.)  phi_cap_s = (1./pi_double)*((sinh(alfa*sqrt(P)))/sqrt(P));  
            if(P<0.)  phi_cap_s = (1./pi_double)*((sin(alfa*sqrt(-P)))/sqrt(-P));   
            if(P==0.) phi_cap_s = alfa/pi_double;

合計時間の約 40% を占めます。

        PP = modulo((r_cc_points + m -K ),(cc*N)); 
        temp.x = temp.x+phi_cap_s*Uj[PP].x; 
        temp.y = temp.y+phi_cap_s*Uj[PP].y;

これは約60％かかります。Visual Profiler により、前者のパフォーマンスがifステートメントの存在に影響されないことを確認しました。倍精度が必要なので、 __exp() ソリューションを避けていることに注意してください。後者の場合、「ランダムな」メモリアクセス Uj[PP] がそれだけの計算パーセンテージを占める可能性があると思います。計算時間を短縮するためのトリックやコメントに関する提案はありますか? 前もって感謝します。

コメントと回答に続くバージョン

回答とコメントで親切に提供された提案に従って、以下のコードになりました。

__global__ void interpolation(cufftDoubleComplex *Uj, double *points, cufftDoubleComplex *result, int N, int M)
{
    int i = threadIdx.x + blockDim.x * blockIdx.x;

    int PP;
    double P,tempd;
    const double alfa=(2.-1./cc)*pi_double-0.01;
    cufftDoubleComplex temp = make_cuDoubleComplex(0.,0.);

    double cc_points=cc*points[i];
    double r_cc_points=rint(cc_points);

    cufftDoubleComplex rtemp[(2*K+1)];
    double phi_cap_s[2*K+1];

    if(i<M) {   
     #pragma unroll //unroll the loop
     for(int m=0; m<(2*K+1); m++) {
         PP = modulo(((int)r_cc_points + m -K ),(cc*N)); 
            rtemp[m] = Uj[PP]; //2

         P = (K*K-(cc_points-(r_cc_points+(double)(m-K)))*(cc_points-(r_cc_points+(double)(m-K))));
         if(P<0.) {tempd=rsqrt(-P); phi_cap_s[m] = (1./pi_double)*((sin(alfa/tempd))*tempd);  }
         else if(P>0.) {tempd=rsqrt(P); phi_cap_s[m] = (1./pi_double)*((sinh(alfa/tempd))*tempd); }
         else phi_cap_s[m] = alfa/pi_double;  
     }

     #pragma unroll //unroll the loop
     for(int m=0; m<(2*K+1); m++) {
         temp.x = temp.x+phi_cap_s[m]*rtemp[m].x; 
           temp.y = temp.y+phi_cap_s[m]*rtemp[m].y; 
     } 

     result[i] = temp; 
     }
 }

特に: 1) グローバル メモリ変数 Uj をサイズ 2*K+1 のレジスタ rtemp 配列に移動しました (私の場合、K は 6 に等しい定数です)。2) 変数 phi_cap_s を 2*K+1 サイズのレジスタに移動しました。3) 以前に使用した 3 つの if ステートメントの代わりに、if ... else ステートメントを使用しました (条件 P<0. と P>0. の発生確率は同じです)。3) 平方根用に追加の変数を定義しました。4) sqrt の代わりに rsqrt を使用しました (私の知る限り、sqrt() は CUDA によって 1/rsqrt() として計算されます)。

それぞれの新機能を 1 つずつ追加し、元のバージョンに対する改善を検証しましたが、どれも関連する改善をもたらさなかったと言わざるを得ません。

実行速度は、1) sin/sinh 関数の計算 (時間の約 40%) によって制限されます。組み込みの数学を「開始推測」として何らかの形で利用することにより、倍精度演算でそれらを計算する方法はありますか? 2) マッピング インデックス PP が原因で、多くのスレッドが同じグローバル メモリ位置 Uj[PP] にアクセスすることになるという事実。これを回避する 1 つの可能性は共有メモリを使用することですが、これは強力なスレッド連携を意味します。

私の質問はです。私は終わりましたか？つまり、コードを改善する手段はありますか? NVIDIA Visual Profiler でコードをプロファイリングした結果は次のとおりです。

IPC = 1.939 (compute capability 2.1);
Global Memory Load Efficiency = 38.9%;
Global Memory Store Efficiency = 18.8%;
Warp Execution Efficiency = 97%;
Instruction Replay Overhead = 0.7%;

最後に、この議論は CUDA での議論にリンクしていることに注意したいと思います: CUDA における 1 次元 3 次スプライン補間

共有メモリを使用するバージョン

共有メモリの使用に関する実現可能性調査を行いました。N=64全体Ujが共有メモリに収まるように配慮しました。以下はコードです（基本的には私のオリジナルバージョンです）

    __global__ void interpolation_shared(cufftDoubleComplex *Uj, double *points, cufftDoubleComplex *result, int N, int M)
 {
         int i = threadIdx.x + blockDim.x * blockIdx.x;

     int PP;
     double P;
     const double alfa=(2.-1./cc)*pi_double-0.01;
     double phi_cap_s;
     cufftDoubleComplex temp;

     double cc_points=cc*points[i];
     double r_cc_points=rint(cc*points[i]);

     temp = make_cuDoubleComplex(0.,0.);

     __shared__ cufftDoubleComplex Uj_shared[128];

     if (threadIdx.x < cc*N) Uj_shared[threadIdx.x]=Uj[threadIdx.x];

     if(i<M) {  
         for(int m=0; m<(2*K+1); m++) {
         P = (K*K-(cc_points-(r_cc_points+m-K))*(cc_points-(r_cc_points+m-K)));

         if(P>0.)  phi_cap_s = (1./pi_double)*((sinh(alfa*sqrt(P)))/sqrt(P));  
         if(P<0.)  phi_cap_s = (1./pi_double)*((sin(alfa*sqrt(-P)))/sqrt(-P));  
         if(P==0.) phi_cap_s = alfa/pi_double;        

         PP = modulo((r_cc_points + m -K ),(cc*N)); 
         temp.x = temp.x+phi_cap_s*Uj_shared[PP].x; 
         temp.y = temp.y+phi_cap_s*Uj_shared[PP].y; 
      } 

      result[i] = temp; 
    }
 }

入力配列のサイズが小さいことにもよるかもしれませんが、ここでも結果は大幅には改善されません。

詳細な PTXAS 出力

ptxas : info : Compiling entry function '_Z13interpolationP7double2PdS0_ii' for 'sm_20'
ptxas : info : Function properties for _Z13interpolationP7double2PdS0_ii
  352 bytes stack frame, 0 bytes spill stores, 0 bytes spill loads
ptxas : info : Used 55 registers, 456 bytes cumulative stack size, 52 bytes cmem[0]

最初のワープと m=0 の場合の P の値

 0.0124300933082964
 0.0127183892149176
 0.0135847002913749
 0.0161796378170038
 0.0155488126345702
 0.0138890822153499
 0.0121163187739057
 0.0119998374528905
 0.0131600831194518
 0.0109574866163769
 0.00962949548477354
 0.00695850974164358
 0.00446426651940612
 0.00423369284281705
 0.00632921297092537
 0.00655137618976198
 0.00810202954519923
 0.00597974034698723
 0.0076811348379735
 0.00604267951733561
 0.00402922460255439
 0.00111841719893846
 -0.00180949615796777
 -0.00246283218698551
 -0.00183256444286428
 -0.000462696661685413
 0.000725108980390132
 -0.00126793006072035
 0.00152263101649197
 0.0022499598348702
 0.00463681632275836
 0.00359856091027666

モジュロ関数

__device__ int modulo(int val, int modulus)
{
   if(val > 0) return val%modulus;
   else
   {
       int P = (-val)%modulus;
       if(P > 0) return modulus -P;
       else return 0;
   }
}

答えに従って最適化された剰余関数

__device__ int modulo(int val, int _mod)
{
    if(val > 0) return val&(_mod-1);
    else
    {
        int P = (-val)&(_mod-1);
        if(P > 0) return _mod -P;
        else return 0;
    }
}
4

2 に答える 2

1 
//your code above
cufftDoubleComplex rtemp[(2*K+1)] //if it fits into available registers, assumes K is a constant

if(i<M) {   
#pragma unroll //unroll the loop
    for(int m=0; m<(2*K+1); m++) {

        PP = modulo((r_cc_points + m -K ),(cc*N)); 
        rtemp[m] = Uj[PP]; //2
    }
#pragma unroll
    for(nt m=0; m<(2*K+1); m++) {
        P = (K*K-(cc_points-(r_cc_points+m-K))*(cc_points-(r_cc_points+m-K)));
        // 1
        if(P>0.)  phi_cap_s = (1./pi_double)*((sinh(alfa*sqrt(P)))/sqrt(P));  
        else if(P<0.)  phi_cap_s = (1./pi_double)*((sin(alfa*sqrt(-P)))/sqrt(-P));   
        else phi_cap_s = alfa/pi_double;  

        temp.x = temp.x+phi_cap_s*rtemp[m].x; //3
        temp.y = temp.y+phi_cap_s*rtemp[m].y; 
    }

    result[i] = temp; 
}

説明

  1. これらの条件は相互に排他的であるため、else if と else を追加しました。可能であれば、発生確率の後にステートメントを並べる必要があります。たとえば、P<0 の場合。ほとんどの場合、最初にそれを評価する必要があります。

  2. これにより、要求されたメモリが複数のレジスタにフェッチされます。以前に行ったことにより、計算に間に合うようにメモリを使用できないため、そのスレッドでブロックが発生した可能性があります。また、ワープで 1 つのスレッドがブロックされると、ワープ全体がブロックされることに注意してください。準備完了キューに十分なワープがない場合、ワープの準備が整うまでプログラムはブロックされます。

  3. 不正なメモリ アクセスを補正するために、計算をさらに前に進めました。うまくいけば、以前に行われた計算によって不正なアクセス パターンが補正されます。

これが機能する理由は次のとおりです。

GMEM にあるメモリからのリクエストは、約 >~400 ～ 600 ティックです。スレッドが、その時点で使用できないメモリに対して操作を実行しようとすると、ブロックされます。つまり、各メモリ要求が L1-L2 に存在しない場合、各ワープは続行できるようになるまでその時間以上待機する必要があります。

私が疑っているtemp.x+phi_cap_s*Uj[PP].xのは、まさにそれをやっているということです。各メモリー転送をレジスターにステージング (ステップ 2) し、次のステージに進むことで、メモリーが転送されている間に他の作業を実行できるようになり、レイテンシーを隠すことができます。

手順 3 に到達するまでに、メモリが使用可能になるか、待機時間が短縮されることが期待されます。

rtemp100% の占有率を達成するためにレジスターに収まらない場合は、バッチで実行する必要がある場合があります。

phi_cap_s次のように、配列にして最初のループに入れることもできます。

#pragma unroll //unroll the loop
    for(int m=0; m<(2*K+1); m++) {
        //stage memory first
        PP = modulo((r_cc_points + m -K ),(cc*N)); 
        rtemp[m] = Uj[PP]; //2

        P = (K*K-(cc_points-(r_cc_points+m-K))*(cc_points-(r_cc_points+m-K)));
        // 1
        if(P>0.)  phi_cap_s[m] = (1./pi_double)*((sinh(alfa*sqrt(P)))/sqrt(P));  
        else if(P<0.)  phi_cap_s[m] = (1./pi_double)*((sin(alfa*sqrt(-P)))/sqrt(-P));   
        else phi_cap_s[m] = alfa/pi_double; 

    }
#pragma unroll
    for(nt m=0; m<(2*K+1); m++) {
        temp.x = temp.x+phi_cap_s[m]*rtemp[m].x; //3
        temp.y = temp.y+phi_cap_s[m]*rtemp[m].y; 
    }

編集

表現

P = (K*K-(cc_points-(r_cc_points+(double)(m-K)))*(cc_points-(r_cc_points+(double)(m-K))));

次のように分類できます。

const double cc_diff = cc_points-r_cc_points;
double exp = cc_diff - (double)(m-K);
exp *= exp;
P = (K*K-exp);

これにより、使用される命令の数が減る可能性があります。

編集 2

__global__ void interpolation(cufftDoubleComplex *Uj, double *points, cufftDoubleComplex *result, int N, int M)
{
    int i = threadIdx.x + blockDim.x * blockIdx.x;

    int PP;
    double P,tempd;


    cufftDoubleComplex rtemp[(2*K+1)];
    double phi_cap_s[2*K+1];

    if(i<M) {
         const double cc_points=cc*points[i];
         cufftDoubleComplex temp = make_cuDoubleComplex(0.,0.);

         const double alfa=(2.-1./cc)*pi_double-0.01;


         const double r_cc_points=rint(cc_points);
         const double cc_diff = cc_points-r_cc_points;

     #pragma unroll //unroll the loop
         for(int m=0; m<(2*K+1); m++) {
             PP = m-k; //reuse PP
             double exp = cc_diff - (double)(PP); //stage exp to be used later, will explain

             PP = modulo(((int)r_cc_points + PP ),(cc*N)); 
             rtemp[m] = Uj[PP]; //2


             exp *= exp;
             P = (K*K-exp);

             if(P<0.) {tempd=rsqrt(-P); phi_cap_s[m] = (1./pi_double)*((sin(alfa/tempd))*tempd);  }
             else if(P>0.) {tempd=rsqrt(P); phi_cap_s[m] = (1./pi_double)*((sinh(alfa/tempd))*tempd); }
             else phi_cap_s[m] = alfa/pi_double;  
         }

     #pragma unroll //unroll the loop
         for(int m=0; m<(2*K+1); m++) {
             temp.x = temp.x+phi_cap_s[m]*rtemp[m].x; 
             temp.y = temp.y+phi_cap_s[m]*rtemp[m].y; 
         } 

     result[i] = temp; 
     }
 }

私が行ったことは、if ステートメント内のすべての計算に移動して、計算とメモリ フェッチの両方の観点からリソースを解放することです。最初の if ステートメントでの分岐がわかりませんif(i<M)。m-Kコードに 2 回出てきたように、最初に とを計算PPするときに使用するために入れました。 expPP

他にできることは、変数を設定する場合、その変数の次の使用の間にできるだけ多くの命令を作成するように命令を試して順序付けることです。登録します。したがって、定数 cc_diff を一番上に置きますが、これは命令の 1 つにすぎないため、何の利点も示さない可能性があります。

モジュロ関数

__device__ modulo(int val, int _mod) {
    int p = (val&(_mod-1));// as modulo is always the power of 2
    if(val < 0) {
        return _mod - p;
    } else {
        return p;
    }
}

_mod常に 2 の累乗の整数 ( ) を使用してcc = 2, N = 64, cc*N = 128いるため、mod 演算子の代わりにこの関数を使用できます。これは「はるかに」高速になるはずです。ただし、算術が正しいのでチェックしてください。これはCuda の最適化 - パート II Nvidiaの 14 ページからのものです。

于 2012-12-27T22:29:11.617 に答える
0

おそらく調べたい最適化の 1 つは、高速計算を使用することです。組み込み数学関数を使用し、-use-fast-math オプションでコンパイルします。

組み込み数学

于 2012-12-26T01:48:28.357 に答える