推力の複素数を使ったcuBLASの使用

Question

私のコードでは推力ライブラリの複素数を持つ配列を使用していますので、配列を転置するためにcublasZgeam（）を使いたいと思います。

cuComplex.hからの複素数の使用は、配列に対して多くの計算を行い、cuComplexは* + =などの演算子を定義していないため、好ましいオプションではありません。コンパイル時に

#include "cuComplex.hpp" 

doesntのは、私が言及した演算子を使用することができます：

この

は私が、私はこの https://github.com/jtravs/cuda_complexを発見した

thrust::complex<float> u[xmax][xmax]; 

を移調したいアレイを定義しますが、そのように使用する方法ですwith nvcc

error: no operator "+=" matches these operands 
     operand types are: cuComplex += cuComplex 

いくつかの方法がありますか？ githubからのコードは古いですし、問題があるかもしれませんが、おそらく私が間違って使っているかもしれません。

EDIT：ここでは動作するコードはありますが、talonmiesコードとの違いは単純なカーネルと同じデータへのポインタですが、あなたの抗議にもかかわらず

#include <iostream> 
#include <thrust/fill.h> 
#include <thrust/complex.h> 
#include <cublas_v2.h> 

using namespace std; 

__global__ void test(thrust::complex<double>* u) { 

    u[0] += thrust::complex<double>(3.3,3.3); 
} 

int main() 
{ 
    int xmax = 100; 
    thrust::complex<double> u[xmax][xmax]; 
    double arrSize = sizeof(thrust::complex<double>) * xmax * xmax; 

    thrust::fill(&u[0][0], &u[0][0] + (xmax * xmax), thrust::complex<double>(1.0,1.0)); 
    u[49][51] += thrust::complex<double>(665.0,665.0); 
    u[51][49] *= 2.0; 

    cout << "Before:" << endl; 
    cout << u[49][51] << endl; 
    cout << u[51][49] << endl; 
    cout << u[0][0] << endl; 

    thrust::complex<double> alpha(1.0, 0.0); 
    thrust::complex<double> beta(0.0, 0.0); 
    cublasHandle_t handle; 
    cublasCreate(&handle); 

    cuDoubleComplex* d_u; 
    cuDoubleComplex* d_v; 
    cuDoubleComplex* _alpha = reinterpret_cast<cuDoubleComplex*>(&alpha); 
    cuDoubleComplex* _beta = reinterpret_cast<cuDoubleComplex*>(&beta); 
    cudaMalloc(&d_u, arrSize); 
    cudaMalloc(&d_v, arrSize); 
    cudaMemcpy(d_u, &u[0][0], arrSize, cudaMemcpyHostToDevice); 
    thrust::complex<double>* d_vTest = reinterpret_cast<thrust::complex<double>* >(d_v); 
    cublasZgeam(handle, CUBLAS_OP_T, CUBLAS_OP_N, xmax, xmax, 
        _alpha, d_u, xmax, 
        _beta, d_u, xmax, 
        d_v, xmax); 
    test<<<1,1>>>(d_vTest); 
    cudaMemcpy(u, d_v, arrSize, cudaMemcpyDeviceToHost); 
    cout << "After:" << endl; 
    cout << u[0][0] << endl; 
    cout << u[49][51] << endl; 
    cout << u[51][49] << endl; 

    return 0; 
} 

Answer 1

反して、C++標準ライブラリcomplex（またはthrust::complex）に最も確かCUBLASで動作ありません。 cuComplexおよびcuDoubleComplexは、標準のホスト複合型とバイナリ互換性があるため、デバイス上で複雑なデータを使用するCUBLAS関数に渡すとデータが変換されないように設計されています。

#include <algorithm> 
#include <iostream> 
#include <complex> 
#include <cublas_v2.h> 

using namespace std; 

int main() 
{ 
    int xmax = 100; 
    complex<double> u[xmax][xmax]; 
    double arrSize = sizeof(complex<double>) * xmax * xmax; 

    fill(&u[0][0], &u[0][0] + (xmax * xmax), complex<double>(1.0,1.0)); 
    u[49][51] += complex<double>(665.0,665.0); 
    u[51][49] *= 2.0; 

    cout << "Before:" << endl; 
    cout << u[49][51] << endl; 
    cout << u[51][49] << endl; 

    complex<double> alpha(1.0, 0.0); 
    complex<double> beta(0.0, 0.0); 
    cublasHandle_t handle; 
    cublasCreate(&handle); 

    cuDoubleComplex* d_u; 
    cuDoubleComplex* d_v; 
    cuDoubleComplex* _alpha = reinterpret_cast<cuDoubleComplex*>(&alpha); 
    cuDoubleComplex* _beta = reinterpret_cast<cuDoubleComplex*>(&beta); 
    cudaMalloc(&d_u, arrSize); 
    cudaMalloc(&d_v, arrSize); 
    cudaMemcpy(d_u, &u[0][0], arrSize, cudaMemcpyHostToDevice); 
    cublasZgeam(handle, CUBLAS_OP_T, CUBLAS_OP_N, xmax, xmax, 
        _alpha, d_u, xmax, 
        _beta, d_u, xmax, 
        d_v, xmax); 

    cudaMemcpy(u, d_v, arrSize, cudaMemcpyDeviceToHost); 

    cout << "After:" << endl; 
    cout << u[49][51] << endl; 
    cout << u[51][49] << endl; 

    return 0; 
} 

が構築されており、そのように実行します：

~/SO$ nvcc -std=c++11 -arch=sm_52 -o complex_transpose complex_transpose.cu -lcublas 
~/SO$ ./complex_transpose 
Before: 
(666,666) 
(2,2) 
After: 
(2,2) 
(666,666) 

必要なだけの変更の明示的なキャストです

あなたはコメントで掲示コードに簡単な変更は、ご想像のとおりに正確に動作しますstd::complex<double>タイプからcuDoubleComplexまでです。それとすべてが期待どおりに機能します。

使用推力は、コードがほぼ同じになります。おそらく、

#include <iostream> 
#include <thrust/fill.h> 
#include <thrust/complex.h> 
#include <cublas_v2.h> 

using namespace std; 

int main() 
{ 
    int xmax = 100; 
    thrust::complex<double> u[xmax][xmax]; 
    double arrSize = sizeof(thrust::complex<double>) * xmax * xmax; 

    thrust::fill(&u[0][0], &u[0][0] + (xmax * xmax), thrust::complex<double>(1.0,1.0)); 
    u[49][51] += thrust::complex<double>(665.0,665.0); 
    u[51][49] *= 2.0; 

    cout << "Before:" << endl; 
    cout << u[49][51] << endl; 
    cout << u[51][49] << endl; 

    thrust::complex<double> alpha(1.0, 0.0); 
    thrust::complex<double> beta(0.0, 0.0); 
    cublasHandle_t handle; 
    cublasCreate(&handle); 

    cuDoubleComplex* d_u; 
    cuDoubleComplex* d_v; 
    cuDoubleComplex* _alpha = reinterpret_cast<cuDoubleComplex*>(&alpha); 
    cuDoubleComplex* _beta = reinterpret_cast<cuDoubleComplex*>(&beta); 
    cudaMalloc(&d_u, arrSize); 
    cudaMalloc(&d_v, arrSize); 
    cudaMemcpy(d_u, &u[0][0], arrSize, cudaMemcpyHostToDevice); 
    cublasZgeam(handle, CUBLAS_OP_T, CUBLAS_OP_N, xmax, xmax, 
        _alpha, d_u, xmax, 
        _beta, d_u, xmax, 
        d_v, xmax); 

    cudaMemcpy(u, d_v, arrSize, cudaMemcpyDeviceToHost); 

    cout << "After:" << endl; 
    cout << u[49][51] << endl; 
    cout << u[51][49] << endl; 

    return 0; 
} 

近いものをご使用の場合には、前CUBLASコールにいくつかの初期化を実行するカーネルで推力装置コンテナーを使用して：

#include <iostream> 
#include <thrust/device_vector.h> 
#include <thrust/complex.h> 
#include <thrust/execution_policy.h> 
#include <thrust/copy.h> 
#include <cublas_v2.h> 

__global__ void setup_kernel(thrust::complex<double>* u, int xmax) 
{ 
    u[51 + 49*xmax] += thrust::complex<double>(665.0,665.0); 
    u[49 + 51*xmax] *= 2.0; 
} 

int main() 
{ 
    int xmax = 100; 

    thrust::complex<double> alpha(1.0, 0.0); 
    thrust::complex<double> beta(0.0, 0.0); 
    cublasHandle_t handle; 
    cublasCreate(&handle); 

    thrust::device_vector<thrust::complex<double>> d_u(xmax * xmax, thrust::complex<double>(1.0,1.0)); 
    thrust::device_vector<thrust::complex<double>> d_v(xmax * xmax, thrust::complex<double>(0.,0.)); 
    setup_kernel<<<1,1>>>(thrust::raw_pointer_cast(d_u.data()), xmax); 

    cuDoubleComplex* _d_u = reinterpret_cast<cuDoubleComplex*>(thrust::raw_pointer_cast(d_u.data())); 
    cuDoubleComplex* _d_v = reinterpret_cast<cuDoubleComplex*>(thrust::raw_pointer_cast(d_v.data())); 
    cuDoubleComplex* _alpha = reinterpret_cast<cuDoubleComplex*>(&alpha); 
    cuDoubleComplex* _beta = reinterpret_cast<cuDoubleComplex*>(&beta); 

    cublasZgeam(handle, CUBLAS_OP_T, CUBLAS_OP_N, xmax, xmax, 
        _alpha, _d_u, xmax, 
        _beta, _d_u, xmax, 
        _d_v, xmax); 

    thrust::complex<double> u[xmax][xmax]; 

    thrust::copy(d_u.begin(), d_u.end(), &u[0][0]); 
    std::cout << "Before:" << std::endl; 
    std::cout << u[49][51] << std::endl; 
    std::cout << u[51][49] << std::endl; 

    thrust::copy(d_v.begin(), d_v.end(), &u[0][0]); 
    std::cout << "After:" << std::endl; 
    std::cout << u[49][51] << std::endl; 
    std::cout << u[51][49] << std::endl; 

    return 0; 

}