はじめに

以前の投稿で、OpenACCとThrustライブラリを組み合わせ、それぞれで確保したメモリをもう一方の機能で計算に利用する方法を紹介しました。
今回は、NVIDIAが提供する他のCUDAライブラリでの利用方法を紹介いたします。

まずOpenACCで確保したメモリについて、CUDA nativeなポインタに変換する方法をおさらいします。
OpenACCでは通常通りメモリを確保して、use_device(x)を使うことで、そのコードブロック内ではOpenACCで確保されたメモリへのポインタが、CUDA nativeなポインタとして振る舞います。

#pragma acc data copy(x[0:n])
    {
#pragma acc host_data use_device(x)
        {
            cuda_native_function(x)
        }
    }

この機能を使えば、基本的にはCUDA nativeなポインタを受け取れるライブラリならどれでも利用可能です。実装が大変なコードから作成していき、最後にOpenACC + managed memory + CUDAライブラリの例を示し、managed memoryの注意点を述べます。
今回はNVIDIA HPC SDK 23.11のnvc++コンパイラを使用し動作確認をしています。

cuFFTを呼び出す (CUDA native版)

例として、FFTのCUDA実装を提供するcuFFTと連携を行います。簡単な64点のFFTのforward/backwardを計算し、backward結果を入力と比較することで計算結果を検証します。
まずはCUDA nativeで実装したコードを以下に示します。

#include <cufft.h>
#include <cuda_runtime_api.h>

#include <iostream>
#include <cmath>
#include <complex>
#include <vector>

#define CUDA_SAFE_CALL(x) { \
    const auto ret = (x); \
    if (ret != cudaSuccess) { \
      std::fprintf(stderr, "[l. %d] CUDA runtime error = %d\n", __LINE__, ret); \
    } \
  }

#define CUFFT_SAFE_CALL(x) { \
    const auto ret = (x); \
    if (ret != CUFFT_SUCCESS) { \
      std::fprintf(stderr, "[l. %d] cuFFT error = %d\n", __LINE__, ret); \
    } \
  }

constexpr int num_points = 64;

int main()
{
  static_assert(sizeof(std::complex<double>) == sizeof(cufftDoubleComplex));

  // create host data
  std::vector<std::complex<double>> input(num_points), forward(num_points), backward(num_points);
  for (int i = 0 ; i < num_points ; ++i)
  {
    const double theta = (double)i / (double)num_points * M_PI;
    const double real = cos(9.0 * theta) + 0.5 * cos(-20.0 * theta);
    const double imag = sin(9.0 * theta) + 0.5 * sin(-20.0 * theta);
    input[i] = std::complex<double>(real, imag);
  }

  // create and copy device data
  cufftDoubleComplex *d_input, *d_forward, *d_backward;
  CUDA_SAFE_CALL(cudaMalloc(reinterpret_cast<void**>(&d_input), sizeof(cufftDoubleComplex) * num_points));
  CUDA_SAFE_CALL(cudaMalloc(reinterpret_cast<void**>(&d_forward), sizeof(cufftDoubleComplex) * num_points));
  CUDA_SAFE_CALL(cudaMalloc(reinterpret_cast<void**>(&d_backward), sizeof(cufftDoubleComplex) * num_points));
  CUDA_SAFE_CALL(cudaMemcpy(input.data(), d_input, sizeof(cufftDoubleComplex) * num_points, cudaMemcpyDefault));

  // create FFT plan
  cufftHandle plan;
  CUFFT_SAFE_CALL(cufftCreate(&plan));
  CUFFT_SAFE_CALL(cufftPlan1d(&plan, num_points, CUFFT_Z2Z, 1));

  // forward FFT
  CUFFT_SAFE_CALL(cufftExecZ2Z(plan, d_input, d_forward, CUFFT_FORWARD));

  // backward FFT
  CUFFT_SAFE_CALL(cufftExecZ2Z(plan, d_forward, d_backward, CUFFT_INVERSE));

  // result copy to host
  CUDA_SAFE_CALL(cudaMemcpy(d_forward, forward.data(), sizeof(cufftDoubleComplex) * num_points, cudaMemcpyDefault));
  CUDA_SAFE_CALL(cudaMemcpy(d_backward, backward.data(), sizeof(cufftDoubleComplex) * num_points, cudaMemcpyDefault));

  // verify
  for (int i = 0 ; i < num_points ; ++i)
  {
    const auto back = backward[i] / static_cast<double>(num_points);
    const double norm = std::abs(input[i] - back);
    if (norm >= 1e-15)
      std::cerr << "error i=" << i << ": " << input[i] << " != " << back << std::endl;
  }

  // cleanup CUDA
  CUFFT_SAFE_CALL(cufftDestroy(plan));
  CUDA_SAFE_CALL(cudaFree(d_backward));
  CUDA_SAFE_CALL(cudaFree(d_forward));
  CUDA_SAFE_CALL(cudaFree(d_input));

  return 0;
}

上記のコードを-cudalib=cufftオプションをつけてコンパイルします。-cudalib=cufftはNVIDIA HPCコンパイラのオプションで、指定したCUDAライブラリおよび必要なライブラリ一式をリンクします。-Lと-lオプションの代わりになります。

$ nvc++ -std=c++17 -cudalib=cufft cufft.cxx
$ ./a.out
(計算結果が問題なければ出力はありません)

OpenACCで確保したメモリを使ってcuFFTを呼び出す

CUDA nativeで実装したコードのメモリ確保部分をOpenACCに置き換えると、以下のような実装になります。
OpenACC pragmaにはstd::vector<T>のような型は渡せないため、一度raw pointerに変換してからacc data句でのメモリ確保を指定する必要があります。

今回はacc enter dataおよびacc exit dataを使用しました。
enter/exitは単独の構文になるので、enterではcreateまたはcopyinを、exitでは対となるようにdeleteまたはcopyoutを指定します。

 int main()
{
  // create host data
  ...

  // get raw pointer
  auto input_ptr = reinterpret_cast<cufftDoubleComplex*>(input.data());
  auto forward_ptr = reinterpret_cast<cufftDoubleComplex*>(forward.data());
  auto backward_ptr = reinterpret_cast<cufftDoubleComplex*>(backward.data());

#pragma acc enter data copyin(input_ptr[0:num_points]) \
                       create(forward_ptr[0:num_points], backward_ptr[0:num_points])

  // create FFT plan
  cufftHandle plan;
  CUFFT_SAFE_CALL(cufftCreate(&plan));
  CUFFT_SAFE_CALL(cufftPlan1d(&plan, num_points, CUFFT_Z2Z, 1));

#pragma acc host_data use_device(input_ptr, forward_ptr, backward_ptr)
  {
    // forward FFT
    CUFFT_SAFE_CALL(cufftExecZ2Z(plan, input_ptr, forward_ptr, CUFFT_FORWARD));

    // backward FFT
    CUFFT_SAFE_CALL(cufftExecZ2Z(plan, forward_ptr, backward_ptr, CUFFT_INVERSE));
  }

#pragma acc exit data copyout(forward_ptr[0:num_points], backward_ptr[0:num_points])

  // verify
  ...

  // cleanup CUDA
}

先程のコンパイルオプションに-accをつけるだけで、上記のコードはコンパイル可能です。

$ nvc++ -std=c++17 -cudalib=cufft -acc cufft_with_oacc.cxx
$ ./a.out
(計算結果が問題なければ出力はありません)

OpenACC + managed memoryを使用する場合

最後に、acc dataを使わずにmanaged memoryで完結させる方法を紹介いたします。
managed memoryはCPUとGPUどちらからもアクセス可能なメモリで、データコピーが必要なタイミングで都度CPU-GPU間で行われます。
またポインタはCPUとGPUどちらからでもアクセスできるアドレスになるため、CUDAカーネルにそのまま指定可能でacc host_data use_deviceも不要になります。

しかし、cuFFTを含むCUDAライブラリの呼び出し直後にCPUで計算結果にアクセスする場合、殆どのケースで計算完了を待つ必要があります。
コード上では、CUDA Streamを使ってcuFFTの計算完了を待機しています。

先程までに紹介した2つの例では、cuFFTの計算結果を検証するためにGPU to CPUのメモリコピー (cudaMemcpyとacc exit data copyout) を明示的に行っており、これがcuFFT (GPU) とCPU間の同期処理として機能していました。
ですがmanaged memoryでは明示的なメモリコピーがないため、cuFFTが計算している最中にCPUが計算結果を格納するメモリにアクセスしてしまいます。
同期処理がないため、CPUから見るとメモリの内容は不定です。
配列をsequentialにアクセスするだけならば、managed memoryにアクセスすることで発生する自動的なメモリコピーを同期処理として使えますが、FFTはそうではないため計算全体が完了するのを待つ必要があります。

#include <cufft.h>
#include <cuda_runtime_api.h>
#include <iostream>
#include <cmath>
#include <complex>
#include <vector>

#define CUDA_SAFE_CALL(x) { \
    const auto ret = (x); \
    if (ret != cudaSuccess) { \
      std::fprintf(stderr, "[l. %d] CUDA runtime error = %d\n", __LINE__, ret); \
    } \
  }

#define CUFFT_SAFE_CALL(x) { \
    const auto ret = (x); \
    if (ret != CUFFT_SUCCESS) { \
      std::fprintf(stderr, "[l. %d] cuFFT error = %d\n", __LINE__, ret); \
    } \
  }

constexpr int num_points = 64;

int main()
{
  static_assert(sizeof(std::complex<double>) == sizeof(cufftDoubleComplex));

  // create host data
  std::vector<std::complex<double>> input(num_points), forward(num_points), backward(num_points);
  for (int i = 0 ; i < num_points ; ++i)
  {
    const double theta = (double)i / (double)num_points * M_PI;
    const double real = cos(9.0 * theta) + 0.5 * cos(-20.0 * theta);
    const double imag = sin(9.0 * theta) + 0.5 * sin(-20.0 * theta);
    input[i] = std::complex<double>(real, imag);
  }

  // create FFT plan
  cufftHandle plan;
  CUFFT_SAFE_CALL(cufftCreate(&plan));
  CUFFT_SAFE_CALL(cufftPlan1d(&plan, num_points, CUFFT_Z2Z, 1));

  // forward FFT
  CUFFT_SAFE_CALL(cufftExecZ2Z(plan, reinterpret_cast<cufftDoubleComplex*>(input.data()),
    reinterpret_cast<cufftDoubleComplex*>(forward.data()), CUFFT_FORWARD));

  // backward FFT
  CUFFT_SAFE_CALL(cufftExecZ2Z(plan, reinterpret_cast<cufftDoubleComplex*>(forward.data()),
    reinterpret_cast<cufftDoubleComplex*>(backward.data()), CUFFT_INVERSE));

  // wait all computations (all kernels run on default stream = `0`)
  CUDA_SAFE_CALL(cudaStreamSynchronize(0));

  // verify
  for (int i = 0 ; i < num_points ; ++i)
  {
    const auto back = backward[i] / static_cast<double>(num_points);
    const double norm = std::abs(input[i] - back);
    if (norm >= 1e-15)
      std::cerr << "error i=" << i << ": " << input[i] << " != " << back << std::endl;
  }

  // cleanup
  CUFFT_SAFE_CALL(cufftDestroy(plan));

  return 0;
}

まとめ

今回はOpenACC + CUDAライブラリの連携について説明し、OpenACCとcuFFTの連携によってコードがどう変わるかと、managed memoryを使ったときの注意点を紹介しました。

複雑なメモリアクセスを要する計算かつCPUとGPUを頻繁に往復するデータの場合、managed memoryによる暗黙的なメモリコピーだけではバグが発生する可能性があるため、cudaStreamSynchronizeやcudaDeviceSynchronizeなどを適切に用いて同期が必要になると予想されますので注意してください。