MAGMAでC++のstd::complex<T>をmagmaFloatComplex/magmaDoubleComplexとして使う #C++

概要

MAGMA(GPU/CUDA版)のmagmaFloatComplexやmagmaDoubleComplex(実態は CUDA CのcuFloatComplexやcuDoubleComplex)はstd::complex<float>やstd::complex<double>に暗黙の変換を行うことはできない．
しかし，両者のバイナリレベルでの互換性はあるので，C++のreinterpret_castを使用することで，互いにポインタを変換して使用する事ができる．
従って，MAGMAの複素行列(複素数ポインタ)を引数に取る関数を，

magmaDoubleComplex *d_a;

と宣言したポインタではなく，

std::complex<double> *d_a;

と宣言したポインタを用いて，

magma_zhoge(…, reinterpret_cast<magmaDoubleComplex*>(d_a), …);

と使うことができる．
ポインターが変換できるので，スカラーの場合も，

std::complex<double> a(1., 2.);
magmaDoubleComplex b = *reinterpret_cast<magmaDoubleComplex*>(&a);

などとして，変換することができる．

検証

実験環境

GCC 4.4.7 (-std=c++0x)
CUDA 5.5.22
Intel MKL 11.1
GPU: NVIDIA M2090
MAGMA 1.5.0 Beta2

実験1

以下では，倍精度複素数型のベクトルをコピーする関数を例に，複素数型の挙動について検証していく．

Host側の変数をstd::complex<double>，Device側の変数をmagmaDoubleComplexで定義する．
Host-Device間の転送は，cuBLASのcublasSetVector/cublasGetVector関数を使用する．

コードは以下の通りである．
(エラーハンドリングは冗長なため，省略してあります)

complex_pointer_test.cpp

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

int main()
{
    std::complex<double> *a, *b; // host
    magmaDoubleComplex *d_a, *d_b; // device
    const size_t n = 10;
    const size_t ldda = ((n+31)/32)*32;

    magma_init();

    // alloc
    a = new std::complex<double>[n];
    b = new std::complex<double>[n];
    cudaMalloc((void**)&d_a, ldda*sizeof(magmaDoubleComplex));
    cudaMalloc((void**)&d_b, ldda*sizeof(magmaDoubleComplex));

    // init a and b
    for (int i=0; i<n; i++)
        a[i] = std::complex<double>(i+1., 0.);
    std::fill(b, b+n, std::complex<double>(0., 0.));

    // output a and b
    std::cout << "a: ";
    for (int i=0; i<n; i++)
        std::cout << a[i] << ' ';
    std::cout << std::endl;
    std::cout << "b: ";
    for (int i=0; i<n; i++)
        std::cout << b[i] << ' ';
    std::cout << std::endl;

    // transfer a to d_a
    cublasSetVector(n, sizeof(std::complex<double>), a, 1, d_a, 1);

    // copy d_a to d_b
    magma_zcopy(n, d_a, 1, d_b, 1);

    // transfer d_b to b
    cublasGetVector(n, sizeof(std::complex<double>), d_b, 1, b, 1);

    // output b
    std::cout << "b: ";
    for (int i=0; i<n; i++)
        std::cout << b[i] << ' ';
    std::cout << std::endl;

    // release
    delete[] a;
    delete[] b;
    cudaFree(d_a);
    cudaFree(d_b);

    magma_finalize();

    return 0;
}

コンパイルして実行すると，以下の出力を得た．

a: (1,0) (2,0) (3,0) (4,0) (5,0) (6,0) (7,0) (8,0) (9,0) (10,0) 
b: (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) 
b: (1,0) (2,0) (3,0) (4,0) (5,0) (6,0) (7,0) (8,0) (9,0) (10,0)

実験2

Host側/Device側共に，変数をstd::complex<double>で定義する．
Host-Device間の転送は，cuBLASのcublasSetVector/cublasGetVector関数を使用する．

コードは以下の通りである．
(説明に必要な部分のみ)

complex_pointer_test2.cpp

// …

int main()
{
    std::complex<double> *a, *b; // host
    std::complex<double> *d_a, *d_b; // device

    // …

    cudaMalloc((void**)&d_a, ldda*sizeof(std::complex<double>));
    cudaMalloc((void**)&d_b, ldda*sizeof(std::complex<double>));

    // …

    // copy d_a to d_b
    magma_zcopy(n, d_a, 1, d_b, 1);

    // …
}

コンパイルしようとすると，以下のエラーが出力された．
(magma_zcopyのところ)

complex_pointer_test2.cpp:40: error: cannot convert ‘std::complex<double>*’ to ‘const magmaDoubleComplex*’ for argument ‘2’ to ‘void magma_zcopy(magma_int_t, const magmaDoubleComplex*, magma_int_t, magmaDoubleComplex*, magma_int_t)’

実験1の結果から，バイナリレベルでの型の互換性はあると考えられるので，

reinterpret_cast<magmaDoubleComplex*>()

を使用して，MAGMAの関数に渡すポインタ変数をキャストする．
修正後コードは以下の通りである．

complex_pointer_test2.cpp

// …

int main()
{
    std::complex<double> *a, *b; // host
    std::complex<double> *d_a, *d_b; // device

    // …

    cudaMalloc((void**)&d_a, ldda*sizeof(std::complex<double>));
    cudaMalloc((void**)&d_b, ldda*sizeof(std::complex<double>));

    // …

    // copy d_a to d_b
    // magma_zcopy(n, d_a, 1, d_b, 1);
    magma_zcopy(n, reinterpret_cast<magmaDoubleComplex*>(d_a), 1, reinterpret_cast<magmaDoubleComplex*>(d_b), 1);

    // …

    return 0;
}

コンパイルして実行すると，以下の結果を得た．

a: (1,0) (2,0) (3,0) (4,0) (5,0) (6,0) (7,0) (8,0) (9,0) (10,0) 
b: (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0) 
b: (1,0) (2,0) (3,0) (4,0) (5,0) (6,0) (7,0) (8,0) (9,0) (10,0)

実験1と同じ結果が得られた．
ここまで，転送のためにcuBLASの関数を使用していたが，このテクニックを使用することで，MAGMAの転送関数(magma_zsetvector/magma_zgetvector, etc)を利用できるようになると考えられる．

実験3

ポインターの変換を行っても正しく動作することが判ったので，ポインタを経由すればスカラーも変換することが可能であると考えられる．

検証用ソースコードは以下の通りである．

complex_scalar_cast.cpp

#include <iostream>
#include <complex>
#include <magma.h>

int main()
{
    std::complex<double> a(1., 2.);
    magmaDoubleComplex b;
    b.x = b.y = 0.;

    std::cout << "b: (" << b.x << ',' << b.y << ')' << std::endl;
    b = *reinterpret_cast<magmaDoubleComplex*>(&a);
    std::cout << "b: (" << b.x << ',' << b.y << ')' << std::endl;

    return 0;
}

コンパイルして実行すると以下の出力が得られた．

b: (0,0)
b: (1,2)

以上より，正しく変換できていることが判った．

参考文献

MAGMA 1.5.0 Beta2 のソースファイル
CUDA 5.5 のヘッダーファイル (cuComplex.h, etc)
MAGMA Forum • View topic - std::complex<double> to cuComplexDouble
(本の虫: C++03とC++11の違い：数値ライブラリー編)