CUDAでホストとデバイスのための効率的な疎テンソルへのハッシュテーブルのアプローチ

Question

CPUやGPUのいずれかで再帰によって生成された値を格納するために使用される多次元疎テンソルで動作する効率的なコードを実装したい。この目標を達成するためには、aligned storageのデータを持つハッシュテーブルが、ストレージとパフォーマンスの間で良好な妥協点を示していると私は推測しています。

私はCPUの実装の最小限のバージョンを持っていますが、コードは以下のとおりです。

私の目的は、GPU上のすべてのカーネルが再帰を満たし、与えられた配列にFの対応する値を選択する位置に格納することです。私はテンソルをハッシュテーブルとして表現することで、最小メモリサイズのデータ​​の表現は、この仮定は正しいですか？

グローバルメモリからデバイスのローカルメモリにデータを転送する際に、高性能とコアレッセンスを実現するためにテンソルをローカルメモリに書き込むにはどうすればよいですか？最終的なアプリケーションのハッシュテーブル用のバッファのサイズは、およそ1から80の間になります。

UPDATE：私は構造H_elementに格納されたキーと値へのアクセスを持っている権利構文式を見つけるために管理することができませんでしたので、私はテンソルのために整列ストレージを使用して問題を経験しました。したがって、私は次の簡略化されたバージョンに切り替えています。

#include <iostream> 
#include <vector> 
#include <random> 
#include <stdio.h> 
#include <algorithm> 
using namespace std; 



struct point{ 
    double x; 
    double y; 
    double z; 

    inline point operator=(double c) { 
    x=c; 
    y=c; 
    z=c; 
    return {x,y,z}; 
    } 

    inline point operator=(point a) { 
    x=a.x; 
    y=a.y; 
    z=a.z; 
    return a; 
    } 

    inline point operator+(point a) { 
    return {a.x+x,a.y+y,a.z+z}; 
    } 

    inline point operator-(point a) { 
    return {a.x-x,a.y-y,a.z-z}; 
    } 

    inline point operator*(double k) { 
    return {k*x,k*y,k*z}; 
    } 



    inline bool operator==(point a) { 
    if (a.x==x && a.y==y && a.z==z) 
     return true; 
    else 
     return false; 
    } 


    double norm2() 
    { 
    return x*x + y*y + z*z; 
    } 

}; 

inline point operator*(double k, point p) { 
    return {k*p.x,k*p.y,k*p.z}; 
} 

inline point operator*(point p,double k) { 
    return {k*p.x,k*p.y,k*p.z}; 
} 



static inline bool is_in(int &na, int &nb, int N) 
{ 
    if(N < 0 || N > (na + nb) || na < 0 || nb < 0) { 
     return false; 
    } else { 
     return true; 
    } 
} 
/// elements of the hash table 
template<typename T> 
struct H_element{ 
    int key; 
    T value; 
}; 

template<typename T> 
struct ascending 
{ 
    inline bool operator() (const H_element<T>& struct1, const H_element<T>& struct2) 
    { 
    return (struct1.key < struct2.key); 
    } 
}; 


template <typename T, size_t dim_na=6,size_t dim_nb=6,size_t dim_N=12> 
struct E_coeff_sparce { 


    enum {LEN1 = dim_na }; 
    enum {LEN2 = dim_nb }; 
    enum {LEN3 = dim_N }; 


    enum {MAX_AM = dim_na-3 }; 



    ///container 
    vector<H_element<T> > data; 


    ///hash function 
    static int Hash_func(int na, int nb, int N) { 
    return (nb*dim_na + na)*dim_nb + N; 
    } 

    ///Map from indices to keys 
    T Binary_Search(int na, int nb, int N) { 
    int key=Hash_func(na,nb,N); 

    int iteration = 0, left = 0, right = data.size()-1, mid; 

    while (left <= right) { 
     iteration++; 
     mid = (int) ((left + right)/2); 
     if (data[mid].key == key) { 
     iteration++; 
    return data[mid].value; 
     } 
     else if (key > data[mid].key) 
    left = mid + 1; 
     else 
    right = mid - 1; 
    } 
    return T(0.0); 

    } 

    T operator()(int na,int nb, int N){ 
    return Binary_Search(na, nb, N); 
    } 


    ///generator of the hash table 
    void Do_recurrence(point A, double alpha_a, int ax, point B, double alpha_b, int bx, bool Laplacian=true) 
    { 
    data.clear(); 

    point R = A - B; 

    int na_max = ax + 1; 
    int nb_max = bx + 1; 

    if(Laplacian){ 
     na_max += 2; 
     nb_max += 2; 
    } 
    int tmax = na_max + nb_max - 1; 

    T a = alpha_a; 
    T b = alpha_b; 
    T p = a + b; 
    T factor = -(a * b/p); 
    // if(is_in(0,0,0)){ 
     int key=Hash_func(0,0,0); 
     H_element<T> E= {key,std::exp(factor*R.x*R.x)}; 
     printf("%d  %e \n",key ,exp(factor*R.x*R.x)); 
     data.push_back(E); 
     //} 
    /// na_=0 
    for(int nb_ = 1; nb_ < nb_max; nb_++){ 
     for(int t = 0; t < tmax; t++){ 

    int na_ = 0; 
    int nb_p = nb_ - 1; 
    int tp = t - 1; 
    int tn = t + 1; 

    double E_na__nb_p_tp = 0.0; 
    if(is_in(na_, nb_p, tp)){ 
     E_na__nb_p_tp = Binary_Search(na_, nb_p, tp); 
    } 

    double E_na__nb_p_t = 0; 
    if(is_in(na_, nb_p, t)) { 
     E_na__nb_p_t = Binary_Search(na_, nb_p, t); 
    } 

    double E_na__nb_p_tn = 0; 
    if(is_in(na_, nb_p, tn)) { 
     E_na__nb_p_tn = Binary_Search(na_, nb_p, tn); 
    } 
    if(is_in(na_,nb_,t)){ 
     H_element<T> E= {Hash_func(na_,nb_,t),1.0/(2*p) * E_na__nb_p_tp + a/p* R.x * E_na__nb_p_t + (t + 1)*E_na__nb_p_tn}; 
     if(E.value != 0.0){ 
     data.push_back(E); 
     std::sort(data.begin(), data.end(), ascending<T>()); 
     //printf("%d  %e \n",E.key ,E.value); 
     } 

    } 
     } 
    } 
    //printf("******************\n"); 
    ///na_>0 
    for(int nb_ = 0; nb_ < nb_max; nb_++) { 
      for(int na_ = 1; na_ < na_max; na_++) { 
       for(int t = 0; t < tmax; t++) { 

        int na_p = na_ - 1; 
        int tp = t - 1; 
        int tn = t + 1; 

        double E_na_p_nb__tp = 0; 
        if(is_in(na_p, nb_, tp)) { 
         E_na_p_nb__tp = Binary_Search(na_p, nb_, tp); 
        } 

        double E_na_p_nb__t = 0; 
        if(is_in(na_p, nb_, t)) { 
         E_na_p_nb__t = Binary_Search(na_p, nb_, t); 
        } 

        double E_na_p_nb__tn = 0; 
        if(is_in(na_p, nb_, tn)) { 
         E_na_p_nb__tn = Binary_Search(na_p, nb_, tn); 
        } 

      if(is_in(na_,nb_,t)){ 
       H_element<T> E= {Hash_func(na_,nb_,t), 1.0/(2*p) * E_na_p_nb__tp - b/p* R.x * E_na_p_nb__t + (t + 1)*E_na_p_nb__tn}; 
       if(E.value != 0.0){ 
      data.push_back(E); 
      std::sort(data.begin(), data.end(), ascending<T>()); 
      //printf("%d  %e \n",E.key ,E.value); 
       } 
      } 
       } 
      } 
     } 
    for(int i=0; i<data.size(); i++) 
     printf("%d  %e \n",data[i].key ,data[i].value); 
    printf("stored: %d \n",data.size()); 
    } 



    friend ostream& operator<<(ostream &strm,E_coeff_sparce<T> ob){ 
    /** Formatted output of the tensor to specified stream 
    */ 
    printf("Unfolded Tensor.\n"); 
    printf("....................\n"); 
    cout<<"na "<<right; 
    cout<<"nb "<<right; 
    cout<<"N  "<<right; 
    cout<<"value "<<endl; 


    int data_PRECISION = 8; 
    int data_WIDTH = 15; 
    strm.setf(ios::showpoint); 
    strm.precision(data_PRECISION); 


    for(int i=0; i< dim_na; i++){ 
     for(int j=0; j< dim_nb; j++){ 
    for(int N=0; N< dim_N; N++){ 
     if(is_in(i,j,N)){ 
     strm<<i<<" "; 

     strm<<j<<" "; 

     strm<<N<<" "; 
     strm<< ob.Hash_func(i, j, N)<<"  "; 
     strm<<right; 
     strm.setf(ios::showpoint); 
     strm.precision(data_PRECISION); 
     strm.width(data_WIDTH); 
     strm<<ob(i,j,N); 
     strm<<endl; 
     } 
    } 
     } 
    } 
    return strm; 
    } 

}; 


template<typename T> 
void f(point A, T a, int na, point B, T b, int nb) 
{ 
    E_coeff_sparce<double> E; 
    E.Do_recurrence(A, a, na, B, b, nb, true); 

    T p = a+b; 
    printf("  RESULTS   \n"); 
    printf("........................\n"); 
    printf("na nb  F \n"); 
    for(int i_na=0; i_na<=na; i_na++) 
    for(int i_nb=0; i_nb<=nb; i_nb++) 
     printf("%d %d %e \n",i_na, i_nb, pow(M_PI/p, 1.5)*E(i_na,i_nb,0)); 

    //return pow(M_PI/p, 1.5)*E(na,nb,0); 
    return; 

} 


int main() { 


    point A = {0.,1.43233673, -0.96104039}; 
    point B = {0.0,0.0,0.24026010}; 

    double alpha_a = 3.425250914; 
    double alpha_b = 5.033151319; 


    E_coeff_sparce<double> E; 
    E.Do_recurrence(A, alpha_a, 2, B, alpha_b, 2, true); 

    cout<<E<<endl; 


    //cout<<"F = "<<<<endl; 
    f(A, alpha_a, 2, B, alpha_b, 2); 
    return 0; 
} 

Answer 1

すべてのコンポーネントが順序付けられるようにあなたのハッシュテーブルアプローチは、実際に疎行列を表すためCOOフォーマットに非常に類似しています。

http://docs.nvidia.com/cuda/cusparse/index.html#coordinate-format-coo

COOフォーマットと比較すると、私はあなたがスパーステンサーを保存するためにCSRやCSC形式を使用することをお勧めします。

http://docs.nvidia.com/cuda/cusparse/index.html#compressed-sparse-row-format-csr

CSRとCSCフォーマットは一般に疎な行列演算に使用されています。テンソルに適用するときは、マトリクス（X dh*w）またはマトリックス（X d*wh）としてテンソル（X w X dh）を記憶することができます。次に、cuSPARSEライブラリによって提供される疎なBLASルーチンを、さらなる操作のために使用することができます。

現時点では、削減のような操作について詳しく説明していないため、実際にはどのフォーマットを使用するのが最終的に決めるのは難しいです。あなたのアプローチはメモリサイズで非常に良いようですが、バイナリ検索はGPUでひどいパフォーマンスを示します。