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


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

```c++
magmaDoubleComplex *d_a;
```

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

```c++
std::complex<double> *d_a;
```

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

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

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

```c++
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の結果から，バイナリレベルでの型の互換性はあると考えられるので，

```c++
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&lt;double&gt; to cuComplexDouble](http://icl.cs.utk.edu/magma/forum/viewtopic.php?f=2&t=902)
- ([本の虫: C++03とC++11の違い：数値ライブラリー編](http://cpplover.blogspot.jp/2014/01/c03c11_22.html))