がGPU上でいくつかの計算でCUDA

Question

とマトリックスの行をスケーリング、Iは行列の行を拡張する必要がある指定された行和のすべての要素に対して1

 
| a1,1 a1,2 ... a1,N | | alpha1*a1,1 alpha1*a1,2 ... alpha1*a1,N | 
| a2,1 a2,2 ... a2,N | => | alpha2*a2,1 alpha2*a2,2 ... alpha2*a2,N | 
| .   . | | .        . | 
| aN,1 aN,2 ... aN,N | | alphaN*aN,1 alphaN*aN,2 ... alphaN*aN,N | 

そう

 
alphai = 1.0/(ai,1 + ai,2 + ... + ai,N) 

私はalphaのベクトルが必要です。スケーリングされた行列と私は可能な限り少数のblas呼び出しでこれをしたいと思います。このコードは、nvidia CUDAハードウェア上で動作する予定です。誰でもこれを行うためのスマートな方法を知っていますか？

Answer 1

単位ベクトルでBLAS gemvを使用すると、必要なスケーリング係数（1/alpha）の逆数のベクトルが得られます。それは簡単な部分です。

標準的なBLASには、使用できるアダマール積和演算子のようなものがないため、要素を横に並べることは少し難しくなります。また、BLASについて言及しているので、私は行列の大規模なオーダー・ストレージを使用していると推測していますが、これは行単位の操作ではそれほど簡単ではありません。 実際には遅いの方法はピッチごとにBLAS scalになりますが、ローごとに1回のBLASコールが必要です。とキャッシュメモリとL1キャッシュのコヒーレンシの影響により、 。

私の提案は、2回目の操作に自分のカーネルを使うことです。ちょうど彼らが一緒に行くようスケーリング、列方向行をステップ実行ブロックの束を持っている

template<typename T> 
__global__ void rowscale(T * X, const int M, const int N, const int LDA, 
          const T * ralpha) 
{ 
    for(int row=threadIdx.x; row<M; row+=gridDim.x) { 
     const T rscale = 1./ralpha[row]; 
     for(int col=blockIdx.x; col<N; col+=blockDim.x) 
      X[row+col*LDA] *= rscale; 
    } 
} 

：それはすべてのこと、複雑な、このようなおそらく唯一のものである必要はありません。任意のサイズの列の主要な順序付けられたマトリックスで動作する必要があります。メモリアクセスは統合する必要がありますが、パフォーマンスの心配度に応じて、試すことができる最適化がいくつかあります。それは少なくとも何をすべきかという一般的な考えを与える。

Answer 2

CUBLAS 5.0（ベクトルで表される）対角行列による行列の乗算であるBLAS状ルーチンと呼ばCUBLAS（タイプ）dgmmを導入しました。

（行を拡大します）左のオプションまたは列をスケールします右のオプションがあります。

詳細については、CUBLAS 5.0ドキュメントを参照してください。

だからあなたの問題では、あなたは左のオプションでGPUと使用cublasdgmm上のすべてのアルファを含むベクトルを作成する必要があります。

Answer 3

私はCUDAスラストのthrust::transformとcuBLASのcublas<t>dgmmの使用を考慮して、上記の回答を例として更新したいと考えています。これはすでに

Reduce matrix rows with CUDA

と

以下

Reduce matrix columns with CUDA

で扱ってきたので、私は、スケーリング係数alphaさんの計算をスキップしています完全な例です：

#include <thrust/device_vector.h> 
#include <thrust/reduce.h> 
#include <thrust/random.h> 
#include <thrust/sort.h> 
#include <thrust/unique.h> 
#include <thrust/equal.h> 

#include <cublas_v2.h> 

#include "Utilities.cuh" 
#include "TimingGPU.cuh" 

/**************************************************************/ 
/* CONVERT LINEAR INDEX TO ROW INDEX - NEEDED FOR APPROACH #1 */ 
/**************************************************************/ 
template <typename T> 
struct linear_index_to_row_index : public thrust::unary_function<T,T> { 

    T Ncols; // --- Number of columns 

    __host__ __device__ linear_index_to_row_index(T Ncols) : Ncols(Ncols) {} 

    __host__ __device__ T operator()(T i) { return i/Ncols; } 
}; 

/***********************/ 
/* RECIPROCAL OPERATOR */ 
/***********************/ 
struct Inv: public thrust::unary_function<float, float> 
{ 
    __host__ __device__ float operator()(float x) 
    { 
     return 1.0f/x; 
    } 
}; 

/********/ 
/* MAIN */ 
/********/ 
int main() 
{ 
    /**************************/ 
    /* SETTING UP THE PROBLEM */ 
    /**************************/ 

    const int Nrows = 10;   // --- Number of rows 
    const int Ncols = 3;   // --- Number of columns 

    // --- Random uniform integer distribution between 0 and 100 
    thrust::default_random_engine rng; 
    thrust::uniform_int_distribution<int> dist1(0, 100); 

    // --- Random uniform integer distribution between 1 and 4 
    thrust::uniform_int_distribution<int> dist2(1, 4); 

    // --- Matrix allocation and initialization 
    thrust::device_vector<float> d_matrix(Nrows * Ncols); 
    for (size_t i = 0; i < d_matrix.size(); i++) d_matrix[i] = (float)dist1(rng); 

    // --- Normalization vector allocation and initialization 
    thrust::device_vector<float> d_normalization(Nrows); 
    for (size_t i = 0; i < d_normalization.size(); i++) d_normalization[i] = (float)dist2(rng); 

    printf("\n\nOriginal matrix\n"); 
    for(int i = 0; i < Nrows; i++) { 
     std::cout << "[ "; 
     for(int j = 0; j < Ncols; j++) 
      std::cout << d_matrix[i * Ncols + j] << " "; 
     std::cout << "]\n"; 
    } 

    printf("\n\nNormlization vector\n"); 
    for(int i = 0; i < Nrows; i++) std::cout << d_normalization[i] << "\n"; 

    TimingGPU timerGPU; 

    /*********************************/ 
    /* ROW NORMALIZATION WITH THRUST */ 
    /*********************************/ 

    thrust::device_vector<float> d_matrix2(d_matrix); 

    timerGPU.StartCounter(); 
    thrust::transform(d_matrix2.begin(), d_matrix2.end(), 
         thrust::make_permutation_iterator(
           d_normalization.begin(), 
           thrust::make_transform_iterator(thrust::make_counting_iterator(0), linear_index_to_row_index<int>(Ncols))), 
         d_matrix2.begin(), 
         thrust::divides<float>()); 
    std::cout << "Timing - Thrust = " << timerGPU.GetCounter() << "\n"; 

    printf("\n\nNormalized matrix - Thrust case\n"); 
    for(int i = 0; i < Nrows; i++) { 
     std::cout << "[ "; 
     for(int j = 0; j < Ncols; j++) 
      std::cout << d_matrix2[i * Ncols + j] << " "; 
     std::cout << "]\n"; 
    } 

    /*********************************/ 
    /* ROW NORMALIZATION WITH CUBLAS */ 
    /*********************************/ 
    d_matrix2 = d_matrix; 

    cublasHandle_t handle; 
    cublasSafeCall(cublasCreate(&handle)); 

    timerGPU.StartCounter(); 
    thrust::transform(d_normalization.begin(), d_normalization.end(), d_normalization.begin(), Inv()); 
    cublasSafeCall(cublasSdgmm(handle, CUBLAS_SIDE_RIGHT, Ncols, Nrows, thrust::raw_pointer_cast(&d_matrix2[0]), Ncols, 
        thrust::raw_pointer_cast(&d_normalization[0]), 1, thrust::raw_pointer_cast(&d_matrix2[0]), Ncols)); 
    std::cout << "Timing - cuBLAS = " << timerGPU.GetCounter() << "\n"; 

    printf("\n\nNormalized matrix - cuBLAS case\n"); 
    for(int i = 0; i < Nrows; i++) { 
     std::cout << "[ "; 
     for(int j = 0; j < Ncols; j++) 
      std::cout << d_matrix2[i * Ncols + j] << " "; 
     std::cout << "]\n"; 
    } 

    return 0; 
} 

Utilities.cuとのファイルがありますhereであり、ここでは省略する。 TimingGPU.cuおよびTimingGPU.cuhは、hereに維持され、同様に省略される。

私はケプラーK20cに上記のコードをテストしており、これらは結果である：私はcublasCreate時間を除いたよ、cuBLASタイミングで

    Thrust  cuBLAS 
2500 x 1250  0.20ms  0.25ms 
5000 x 2500  0.77ms  0.83ms 

。これでも、CUDA Thrustバージョンがより便利になるようです。