2

C ++ AMPを使用した並列計算用のアルゴリズム(格子ボルツマン)を最適化しようとしています。そして、メモリレイアウトを最適化するためのいくつかの提案を探して、構造から別のベクトル(ブロックされたベクトル)に1つのパラメーターを削除すると、約10%の増加が得られることがわかりました。

誰かがこれをさらに改善することができるヒント、または私が考慮すべき何かを手に入れましたか?以下は、タイムステップごとに実行される最も時間のかかる関数と、レイアウトに使用される構造です。

struct grid_cell {
//  int blocked;    // Define if blocked
    float n;        // North
    float ne;       // North-East
    float e;        // East
    float se;       // South-East
    float s;
    float sw;
    float w;
    float nw;
    float c;        // Center
};

int collision(const struct st_parameters param, vector<struct grid_cell> &node, vector<struct grid_cell> &tmp_node, vector<int> &obstacle) {
    int x,y;
    int i = 0;
    float c_sq = 1.0f/3.0f;     // Square of speed of sound
    float w0 = 4.0f/9.0f;       // Weighting factors
    float w1 = 1.0f/9.0f;
    float w2 = 1.0f/36.0f;

    int chunk = param.ny/20;
    float total_density = 0;

    float u_x,u_y;              // Avrage velocities in x and y direction
    float u[9];                 // Directional velocities
    float d_equ[9];             // Equalibrium densities
    float u_sq;                 // Squared velocity
    float local_density;        // Sum of densities in a particular node

    for(y=0;y<param.ny;y++) {
        for(x=0;x<param.nx;x++) {
            i = y*param.nx + x; // Node index
            // Dont consider blocked cells
            if (obstacle[i] == 0) {
                // Calculate local density
                local_density = 0.0;
                local_density += tmp_node[i].n;
                local_density += tmp_node[i].e;
                local_density += tmp_node[i].s;
                local_density += tmp_node[i].w;
                local_density += tmp_node[i].ne;
                local_density += tmp_node[i].se;
                local_density += tmp_node[i].sw;
                local_density += tmp_node[i].nw;
                local_density += tmp_node[i].c;

                // Calculate x velocity component
                u_x = (tmp_node[i].e + tmp_node[i].ne + tmp_node[i].se - 
                      (tmp_node[i].w + tmp_node[i].nw + tmp_node[i].sw)) 
                      / local_density;
                // Calculate y velocity component
                u_y = (tmp_node[i].n + tmp_node[i].ne + tmp_node[i].nw - 
                      (tmp_node[i].s + tmp_node[i].sw + tmp_node[i].se)) 
                      / local_density;
                // Velocity squared
                u_sq = u_x*u_x +u_y*u_y;

                // Directional velocity components;
                u[1] =  u_x;        // East
                u[2] =        u_y;  // North
                u[3] = -u_x;        // West
                u[4] =      - u_y;  // South
                u[5] =  u_x + u_y;  // North-East
                u[6] = -u_x + u_y;  // North-West
                u[7] = -u_x - u_y;  // South-West
                u[8] =  u_x - u_y;  // South-East

                // Equalibrium densities
                // Zero velocity density: weight w0
                d_equ[0] = w0 * local_density * (1.0f - u_sq / (2.0f * c_sq));
                // Axis speeds: weight w1
                d_equ[1] = w1 * local_density * (1.0f + u[1] / c_sq
                                 + (u[1] * u[1]) / (2.0f * c_sq * c_sq)
                                 - u_sq / (2.0f * c_sq));
                d_equ[2] = w1 * local_density * (1.0f + u[2] / c_sq
                                 + (u[2] * u[2]) / (2.0f * c_sq * c_sq)
                                 - u_sq / (2.0f * c_sq));
                d_equ[3] = w1 * local_density * (1.0f + u[3] / c_sq
                                 + (u[3] * u[3]) / (2.0f * c_sq * c_sq)
                                 - u_sq / (2.0f * c_sq));
                d_equ[4] = w1 * local_density * (1.0f + u[4] / c_sq
                                 + (u[4] * u[4]) / (2.0f * c_sq * c_sq)
                                 - u_sq / (2.0f * c_sq));
                // Diagonal speeds: weight w2
                d_equ[5] = w2 * local_density * (1.0f + u[5] / c_sq
                                 + (u[5] * u[5]) / (2.0f * c_sq * c_sq)
                                 - u_sq / (2.0f * c_sq));
                d_equ[6] = w2 * local_density * (1.0f + u[6] / c_sq
                                 + (u[6] * u[6]) / (2.0f * c_sq * c_sq)
                                 - u_sq / (2.0f * c_sq));
                d_equ[7] = w2 * local_density * (1.0f + u[7] / c_sq
                                 + (u[7] * u[7]) / (2.0f * c_sq * c_sq)
                                 - u_sq / (2.0f * c_sq));
                d_equ[8] = w2 * local_density * (1.0f + u[8] / c_sq
                                 + (u[8] * u[8]) / (2.0f * c_sq * c_sq)
                                 - u_sq / (2.0f * c_sq));

                // Relaxation step
                node[i].c = (tmp_node[i].c + param.omega * (d_equ[0] - tmp_node[i].c));
                node[i].e = (tmp_node[i].e + param.omega * (d_equ[1] - tmp_node[i].e));
                node[i].n = (tmp_node[i].n + param.omega * (d_equ[2] - tmp_node[i].n));
                node[i].w = (tmp_node[i].w + param.omega * (d_equ[3] - tmp_node[i].w));
                node[i].s = (tmp_node[i].s + param.omega * (d_equ[4] - tmp_node[i].s));
                node[i].ne = (tmp_node[i].ne + param.omega * (d_equ[5] - tmp_node[i].ne));
                node[i].nw = (tmp_node[i].nw + param.omega * (d_equ[6] - tmp_node[i].nw));
                node[i].sw = (tmp_node[i].sw + param.omega * (d_equ[7] - tmp_node[i].sw));
                node[i].se = (tmp_node[i].se + param.omega * (d_equ[8] - tmp_node[i].se));
            }
        }
    }
    return 1;
}
4

4 に答える 4

7

一般に、異なるCPUで使用されるデータが共有されておらず(簡単)、同じキャッシュライン上にないことを確認する必要があります(偽共有。たとえば、ここを参照してください:偽共有は面白くない)。同じCPUで使用されるデータは、キャッシュの恩恵を受けるために互いに接近している必要があります。

于 2012-03-13T13:24:51.400 に答える
6

現在のGPUは、メモリレイアウトに依存していることで有名です。アプリケーションの詳細がない場合は、次のことを検討することをお勧めします。

  1. ユニットストライドアクセスは非常に重要であるため、GPUは「構造体の配列」よりも「配列の構造体」を優先します。「ブロックされた」フィールドをベクトル「障害物」に移動したように、「grid_cell」のすべてのフィールドを個別のベクトルに変換すると有利なはずです。これは、すべてのフィールドにアクセスしないループの場合にも、CPUにメリットがあることを示しています。

  2. 「障害物」が非常にまばらである場合(これはありそうもないと思います)、それを独自のベクトルに移動することは特に価値があります。CPUなどのGPUは、キャッシュラインまたはその他の形式でメモリシステムから複数のワードをロードし、データの一部が必要ない場合は帯域幅を浪費します。多くのシステムでは、メモリ帯域幅がボトルネックリソースであるため、帯域幅を削減する方法が役立ちます。

  3. これはより推測的ですが、すべての出力ベクトルを書き込んでいるので、メモリサブシステムが単に上書きされる「ノード」の値の読み取りを回避している可能性があります

  4. 一部のシステムでは、オンチップメモリ​​がバンクに分割されており、構造内に奇数のフィールドがあると、バンクの競合を取り除くのに役立つ場合があります。

  5. 一部のシステムはロードとストアも「ベクトル化」するため、構造から「ブロック」を再度削除すると、より多くのベクトル化が可能になる場合があります。配列構造体への移行により、この心配が軽減されます。

C++AMPに関心をお寄せいただきありがとうございます。

デビッドキャラハン

http://blogs.msdn.com/b/nativeconcurrency/ C++AMPチームブログ

于 2012-03-16T18:08:21.953 に答える
1

いくつかの小さな一般的なトップス:

  • 複数のプロセッサ間で共有されるデータ構造は、読み取り専用である必要があります。

  • 変更が必要なデータ構造はプロセッサに固有であり、別のプロセッサが必要とするデータとメモリの局所性を共有しません。

  • コードが連続してスキャンされるようにメモリが配置されていることを確認してください(大きなステップを踏んだり、飛び回ったりしないでください)。

于 2012-03-13T15:47:22.613 に答える
0

このトピックを調べている人のために、いくつかのヒントがあります。Lattice-Boltzmannは、一般的に帯域幅が制限されています。つまり、そのパフォーマンスは、主にメモリからロードおよびメモリに書き込むことができるデータの量に依存します。

  • 非常に効率的なコンパイル型プログラミング言語を使用します。CPUベースの実装には、CまたはC++が適しています。

  • 高帯域幅のアーキテクチャを選択してください。CPUの場合、これは高クロックRAMと多数のメモリチャネル(クアッドチャネル以上)を意味します。

  • これにより、キャッシュのプリフェッチを効果的に利用する適切な線形メモリレイアウトを使用することが重要になります。データは、メモリ内で小さな部分、いわゆるキャッシュラインに配置されます。プロセッサが要素にアクセスするときはいつでも、それが存在するキャッシュライン全体(最新のアーキテクチャでは64バイト)がロードされます。これは、8つのdouble値または16のfloat値が一度にロードされることを意味します。マルチコアプロセッサはL3キャッシュを共有しているため、これが問題になることはありませんが、同じキャッシュラインへの変更を同期させる必要があり、他のプロセッサで問題が発生するため、メニーコアアーキテクチャで問題が発生する可能性があります。別のプロセッサが処理しているデータを処理しています(偽共有)。これは、パディングを導入することで回避できます、つまり、残りのキャッシュ行を埋めるために使用しない要素を追加します。D3Q27-latticeの27の速度で離散化してセルを更新するとし、double(8バイト)の場合、データは4つの異なるキャッシュラインにあります。32バイト(4 * 8バイト)に一致するように、5つのdoubleのパディングを追加する必要があります。

unsigned int const PAD          = (64 - sizeof(double)*D3Q27.SPEEDS % 64); ///< padding: number of doubles
size_t       const MEM_SIZE_POP = sizeof(double)*NZ*NY*NX*(SPEEDS+PAD);    ///< amount of memory to be allocated

ほとんどのコンパイラは、配列の先頭をキャッシュラインに自然に揃えるので、それを処理する必要はありません。

  • 線形インデックスはアクセスに不便です。したがって、アクセスを可能な限り効率的に設計する必要があります。ラッパークラスを作成できます。いずれの場合も、これらの関数をインライン化します。つまり、すべての呼び出しがコード内の定義に置き換えられます。
inline size_t const D3Q27_PopIndex(unsigned int const x, unsigned int const y, unsigned int const z, unsigned int const d)
{
    return (D3Q27.SPEEDS + D3Q27.PAD)*(NX*(NY*z + y) + x) + D3Q27.SPEEDS*p + d;
}
  • さらに、たとえば3次元空間ループブロッキングOpenMPのスケーリングの問題)を使用して、計算と通信の比率を最大化することにより、キャッシュの局所性を高めることができます。つまり、すべてのコードが単一のセルではなくセルのキューブで機能します。

  • 一般に、実装は2つの異なる母集団AとBを利用し、ある実装から別の実装への衝突とストリーミングを実行します。これは、メモリ内のすべての値が、衝突前と衝突後の2回存在することを意味します。ステップを再結合し、1つのポピュレーションコピーをメモリに保持するだけでよいように保存するためのさまざまな戦略があります。たとえば、P。Baileyetal。によって提案されたAAパターンを参照してください。-「グラフィックプロセッサを使用した格子ボルツマン流体シミュレーションの加速」(https://www2.cs.arizona.edu/people/pbailey/Accelerating_GPU_LBM.pdf)または難解なツイストM. Geier&M.Schönherr著-「EsotericTwist:A Efficient in-place Streaming Algorithmus for the Lattice Boltzmann Method on Massively Parallel Hardware」(https://pdfs.semanticscholar.org/ea64/3d63667900b60e6ff49f2746211700e63802.pdf)。私はマクロを使用して最初の実装を行いました。これは、母集団へのすべてのアクセスが次の形式のマクロを呼び出すことを意味します。

#define O_E(a,b)       a*odd + b*(!odd)

#define  READ_f_0      D3Q27_PopIndex(x,           y, z,           0, p)
#define  READ_f_1      D3Q27_PopIndex(O_E(x_m, x), y, z, O_E( 1,  2), p)
#define  READ_f_2      D3Q27_PopIndex(O_E(x_p, x), y, z, O_E( 2,  1), p)
...

#define  WRITE_f_0     D3Q27_PopIndex(x,           y, z,           0, p)
#define  WRITE_f_1     D3Q27_PopIndex(O_E(x_p, x), y, z, O_E( 2,  1), p)
#define  WRITE_f_2     D3Q27_PopIndex(O_E(x_m, x), y, z, O_E( 1,  2), p)
...
  • 複数の相互作用の母集団がある場合は、グリッドのマージを使用します。インデックスをメモリに直線的に配置し、2つの異なる母集団を並べて配置します。母集団pへのアクセスは、次のように機能します。
inline size_t const D3Q27_PopIndex(unsigned int const x, unsigned int const y, unsigned int const z, unsigned int const d, unsigned int const p = 0)
{
    return (D3Q27.SPEEDS*D3Q27.NPOP + D3Q27.PAD)*(NX*(NY*z + y) + x) + D3Q27.SPEEDS*p + d;
}
  • 通常のグリッドの場合、アルゴリズムを可能な限り予測可能にしますすべてのセルに衝突とストリーミングを実行させ、その後、境界条件を逆に実行させます。アルゴリズムに直接寄与しないセルが多数ある場合は、パディングにも保存できる論理マスクを使用してそれらを省略してください。

  • コンパイル時にコンパイラにすべてを知らせます。たとえば、すべての境界条件を書き直す必要がないように、インデックスの変更を処理する関数を使用した境界条件のテンプレート。

  • 最新のアーキテクチャにはSIMD操作を実行できるレジスタがあるため、複数のデータに対して同じ命令を実行します。一部のプロセッサ(AVX-512)は、最大512ビットのデータを処理できるため、32は単一の数値のほぼ2倍の速度です。これは、ギャザリングとスキャッタリングが導入されて以来(https://en.wikipedia.org/wiki/Gather-scatter_(vector_addressing))、特にLBMにとって非常に魅力的なようですが、現在の帯域幅の制限があります(おそらくそれだけの価値があります)。 DDR5と少数のコアを備えたプロセッサでは、これは面倒な価値がないと私は考えています。シングルコアのパフォーマンスと並列スケーリングの方が優れています(M. Wittmann etal。-「パフォーマンス分析のテストベッドとしての格子ボルツマンベンチマークカーネル」-https: //arxiv.org/abs/1711.11468)ただし、帯域幅が制限されているため、アルゴリズム全体のパフォーマンスは向上しません。したがって、帯域幅ではなくコンピューティング能力によって制限されるアーキテクチャでのみ意味があります。Xeon Phiアーキテクチャでは、結果は注目に値するようです。-「XeonPhiプロセッサでの格子ボルツマン法の高性能SIMD実装」(https://onlinelibrary.wiley.com/doi/abs/10.1002/cpe.5072)。

私の意見では、これのほとんどは、単純な2D実装のために努力する価値はありません。そこで簡単な最適化、ループブロッキング、線形メモリレイアウトを実行しますが、より複雑なアクセスパターンは忘れてください。3Dでは、その効果は計り知れません。最新の12コアプロセッサのD3Q19で、最大95%の並列スケーラビリティと150Mlupsを超える全体的なパフォーマンスを達成しました。パフォーマンスを向上させるには、帯域幅に最適化されたCUDACを備えたGPUなどのより適切なアーキテクチャに切り替えます。

于 2019-10-20T12:19:50.053 に答える