読者です 読者をやめる 読者になる 読者になる

koturnの日記

転職したい社会人1年生の技術系日記

OpenCLのオフラインコンパイル

C++ OpenCL

はじめに

OpenCLといえば,カーネルのコードに以下の2つのコンパイル方式がある.

オンラインコンパイルは,実行時にOpenCLカーネルコードを文字列として関数に渡し,プログラムオブジェクトを構築する手法である. 反対に,オフラインコンパイルは事前にOpenCLカーネルコードをコンパイルし,コンパイル結果のバイナリを生成しておく. そして,生成したバイナリを実行時に読み込んで,プログラムオブジェクトを構築する手法である. オフラインコンパイルは事前にコンパイルを行う分,実行時のコンパイル時間を削減することができるわけだ.

オフラインコンパイル

オフラインコンパイルには,大別して2つの方針がある.

  1. OpenCLSDK付属のオフラインコンパイラを利用する
  2. OpenCLAPIを用いて,コンパイルプログラムを書く

前者はSDKが必須となるが,後者は不要となる. この記事では,OpenCLAPIを用いて,オフラインコンパイラを書くことにする.

オフラインコンパイルを行うプログラムを作成する

clCreateProgramWithSource()clBuildProgram() により生成した cl_program から clGetProgramInfo() を用いることで,バイナリの情報およびバイナリそのものを取得する形になる. 具体的には以下のようになる. このコードは,0番目のプラットフォームID,0番目のデバイスIDを対象に,コマンドライン引数で指定したカーネルソースコードのオフラインコンパイルを行うものだ.

  • oclc.cpp
// oclc.cpp
// g++ -gnu++11 -O3 oclc.cpp -lOpenCL -o oclc
#include <iostream>
#include <cmath>
#include <cstdlib>
#include <cstring>
#include <ctime>
#include <fstream>
#include <memory>
#include <vector>

#ifdef __APPLE__
#  include <OpenCL/opencl.h>
#else
#  include <CL/cl.h>
#endif


static constexpr cl_uint kNDefaultPlatformEntry = 16;
static constexpr cl_uint kNDefaultDeviceEntry = 16;


/*!
 * @brief プラットフォームIDを取得
 * @param [in] nPlatformEntry  取得するプラットフォームID数の上限
 * @return  プラットフォームIDを格納した std::vector
 */
static inline std::vector<cl_platform_id>
getPlatformIds(cl_uint nPlatformEntry = kNDefaultPlatformEntry)
{
  std::vector<cl_platform_id> platformIds(nPlatformEntry);
  cl_uint nPlatform;
  if (clGetPlatformIDs(nPlatformEntry, platformIds.data(), &nPlatform) != CL_SUCCESS) {
    std::cerr << "clGetPlatformIDs() failed" << std::endl;
    std::exit(EXIT_FAILURE);
  }
  platformIds.resize(nPlatform);
  return platformIds;
}


/*!
 * @brief デバイスIDを取得
 * @param [in] platformId    デバイスIDの取得元のプラットフォームのID
 * @param [in] nDeviceEntry  取得するデバイスID数の上限
 * @param [in] deviceType    取得対象とするデバイス
 * @return デバイスIDを格納した std::vector
 */
static inline std::vector<cl_device_id>
getDeviceIds(const cl_platform_id& platformId, cl_uint nDeviceEntry = kNDefaultDeviceEntry, cl_int deviceType = CL_DEVICE_TYPE_DEFAULT)
{
  std::vector<cl_device_id> deviceIds(nDeviceEntry);
  cl_uint nDevice;
  if (clGetDeviceIDs(platformId, deviceType, nDeviceEntry, deviceIds.data(), &nDevice) != CL_SUCCESS) {
    std::cerr << "clGetDeviceIDs() failed" << std::endl;
    std::exit(EXIT_FAILURE);
  }
  deviceIds.resize(nDevice);
  return deviceIds;
}


/*!
 * @brief 指定されたファイル名のファイルを読み込み,std::stringに格納して返却する
 * @param [in] filename  読み込むファイル名
 * @return  ファイルの内容を格納したstd::string
 */
static inline std::string
readSource(const std::string& filename) noexcept
{
  std::ifstream ifs(filename.c_str());
  if (!ifs.is_open()) {
    std::cerr << "Failed to open " << filename << std::endl;
    std::exit(EXIT_FAILURE);
  }
  return std::string((std::istreambuf_iterator<char>(ifs)), std::istreambuf_iterator<char>());
}


/*!
 * @brief 指定された複数のファイル名のファイルを読み込み,std::stringに格納して返却する
 * @param [in] filenames  読み込み対象のファイル名を格納したstd::vector
 * @return  ファイルの内容を格納したstd:stringを格納したstd::vector
 */
static inline std::vector<std::string>
readSource(const std::vector<std::string>& filenames)
{
  std::vector<std::string> srcs(filenames.size());
  for (decltype(srcs)::size_type i = 0; i < srcs.size(); i++) {
    srcs[i] = readSource(filenames[i]);
  }
  return srcs;
}


/*!
 * @brief 指定されたファイル名から拡張子を除いた文字列を返却する
 * @param [in] filename  拡張子を取り除きたいファイル名
 * @return  拡張子を除いたファイル名
 */
static inline std::string
removeSuffix(const std::string& filename) noexcept
{
  return filename.substr(0, filename.find_last_of("."));
}


/*!
 * @brief このプログラムのエントリポイント
 * @param [in] argc  コマンドライン引数の数
 * @param [in] argv  コマンドライン引数
 * @return 終了ステータス
 */
int
main(int argc, char* argv[])
{
  // OpenCLのコンパイラに渡すオプション文字列
  // 今回はとりあえず空
  static const char kOptStr[] = "";

  std::vector<std::string> args(argc - 1);
  if (args.size() < 1) {
    std::cerr << "Please specify only one or more source file" << std::endl;
    return EXIT_FAILURE;
  }
  for (decltype(args)::size_type i = 0; i < args.size(); i++) {
    args[i] = std::string(argv[i + 1]);
  }

  // プラットフォームを取得
  std::vector<cl_platform_id> platformIds = getPlatformIds(1);

  // デバイスを取得
  std::vector<cl_device_id> deviceIds = getDeviceIds(platformIds[0], 1, CL_DEVICE_TYPE_DEFAULT);

  // コンテキスト生成
  cl_int errCode;
  std::unique_ptr<std::remove_pointer<cl_context>::type, decltype(&clReleaseContext)> context(
      clCreateContext(nullptr, 1, &deviceIds[0], nullptr, nullptr, &errCode), clReleaseContext);
  if (errCode != CL_SUCCESS) {
    std::cerr << "clCreateContext() failed" << std::endl;
    return EXIT_FAILURE;
  }

  // ソースコード読み込み
  std::vector<std::string> kernelSources = readSource(args);
  std::pair<std::vector<const char*>, std::vector<std::string::size_type> > kernelSourcePairs;
  kernelSourcePairs.first.reserve(kernelSources.size());
  kernelSourcePairs.second.reserve(kernelSources.size());
  for (const auto& kernelSource : kernelSources) {
    kernelSourcePairs.first.emplace_back(kernelSource.c_str());
    kernelSourcePairs.second.emplace_back(kernelSource.length());
  }

  // プログラム生成
  // 複数ソースファイルに対応
  std::unique_ptr<std::remove_pointer<cl_program>::type, decltype(&clReleaseProgram)> program(
      clCreateProgramWithSource(
        context.get(),
        static_cast<cl_uint>(kernelSourcePairs.first.size()),
        kernelSourcePairs.first.data(),
        kernelSourcePairs.second.data(),
        &errCode),
      clReleaseProgram);
  if (errCode != CL_SUCCESS) {
    std::cerr << "clCreateProgramWithSource() failed" << std::endl;
    return EXIT_FAILURE;
  }

  // カーネルソースコードのコンパイル
  switch (clBuildProgram(program.get(), 1, &deviceIds[0], kOptStr, nullptr, nullptr)) {
    case CL_SUCCESS:
      break;
    case CL_BUILD_PROGRAM_FAILURE:
      {
        // コンパイルエラーを表示
        std::array<char, 2048> buildLog;
        std::size_t logSize;
        clGetProgramBuildInfo(program.get(), deviceIds[0], CL_PROGRAM_BUILD_LOG, buildLog.size(), buildLog.data(), &logSize);
        std::cerr << "Compile error:\n" << buildLog.data() << std::endl;
      }
      break;
    case CL_INVALID_BUILD_OPTIONS:
      std::cerr << "Invalid option is specified" << std::endl;
      return EXIT_FAILURE;
    default:
      std::cerr << "clBuildProgram() failed" << std::endl;
      return EXIT_FAILURE;
  }

  // デバイス数を取得 (このプログラムでは1が返却されるはず)
  cl_uint nDevice;
  if (clGetProgramInfo(program.get(), CL_PROGRAM_NUM_DEVICES, sizeof(nDevice), &nDevice, nullptr) !=  CL_SUCCESS) {
    std::cerr << "clGetProgramInfo() failed" << std::endl;
  }

  // 各デバイス向けのコンパイル後のバイナリのサイズを取得
  std::unique_ptr<std::size_t[]> binSizes(new std::size_t[nDevice]);
  if (clGetProgramInfo(program.get(), CL_PROGRAM_BINARY_SIZES, sizeof(std::size_t) * nDevice, binSizes.get(), nullptr) != CL_SUCCESS) {
    std::cerr << "clGetProgramInfo() failed" << std::endl;
  }

  // コンパイル後のバイナリをコピー
  std::vector<std::unique_ptr<char> > bins(nDevice);
  for (std::size_t i = 0; i < nDevice; i++) {
    bins[i] = std::unique_ptr<char>(binSizes[i] == 0 ? nullptr : new char[binSizes[i]]);
  }
  if (clGetProgramInfo(program.get(), CL_PROGRAM_BINARIES, sizeof(char*) * nDevice, bins.data(), nullptr) != CL_SUCCESS) {
    std::cerr << "clGetProgramInfo() failed" << std::endl;
    return EXIT_FAILURE;
  }

  // コピーしたバイナリを全てファイルに出力
  std::string basename = removeSuffix(args[0]);
  for (std::size_t i = 0; i < nDevice; i++) {
    if (bins[i] == nullptr) {
      continue;
    }
    std::string filename = basename + ".bin";
    if (nDevice > 1) {
      filename += "." + std::to_string(i);
    }
    std::ofstream ofs(filename, std::ios::binary);
    if (ofs.is_open()) {
      ofs.write(bins[i].get(), binSizes[i]);
    } else {
      std::cerr << "Failed to open: " << filename << std::endl;
    }
  }

  return EXIT_SUCCESS;
}

コンパイルは以下のようにして行う.

$ g++ -gnu++11 -O3 oclc.cpp -lOpenCL -o oclc

次に,以下のようなカーネルのコードを用意する.

  • kernel.cl
// kernel.cl
__kernel void
vecAdd(__global float* z, __global const float* x, __global const float* y, int n)
{
  const int para = 4;
  const int end = (n / para) * para;

  for (int i = 0; i < end; i += para) {
    float4 vtmp = vload4(0, x + i) + vload4(0, y + i);
    vstore4(vtmp, 0, z + i);
  }

  for (int i = end; i < n; i++) {
    z[i] = x[i] + y[i];
  }
}

そして,以下のようにしてカーネルソースコードコンパイルする.

$ ./oclc kernel.cl

これでカレントディレクトリに kernel.bin が生成されていれば,オフラインコンパイルは成功である. コンパイルエラーがある場合は,エラーメッセージを表示するようにしてある.

生成したカーネルバイナリのテストには以下のコードを用いる.

// main.cpp
// g++ -std=gnu++11 -O3 main.cpp -lOpenCL -o main
#include <cstdlib>
#include <fstream>
#include <iostream>
#include <memory>
#include <random>
#include <string>

#ifdef __APPLE__
#  include <OpenCL/opencl.h>
#else
#  include <CL/cl.h>
#endif


static constexpr cl_uint kNDefaultPlatformEntry = 16;
static constexpr cl_uint kNDefaultDeviceEntry = 16;




/*!
 * @brief アラインメントされたメモリを動的確保する関数
 * @param [in] size       確保するメモリサイズ (単位はbyte)
 * @param [in] alignment  アラインメント (2のべき乗を指定すること)
 * @return  アラインメントし,動的確保されたメモリ領域へのポインタ
 */
template<typename T = void*, typename std::enable_if<std::is_pointer<T>::value, std::nullptr_t>::type = nullptr>
static inline T
alignedMalloc(std::size_t size, std::size_t alignment) noexcept
{
#if defined(_MSC_VER) || defined(__MINGW32__)
  return reinterpret_cast<T>(_aligned_malloc(size, alignment));
#else
  void* p;
  return reinterpret_cast<T>(posix_memalign(&p, alignment, size) == 0 ? p : nullptr);
#endif  // defined(_MSC_VER) || defined(__MINGW32__)
}


/*!
 * @brief アラインメントされたメモリを解放する関数
 * @param [in] ptr  解放対象のメモリの先頭番地を指すポインタ
 */
static inline void
alignedFree(void* ptr) noexcept
{
#if defined(_MSC_VER) || defined(__MINGW32__)
  _aligned_free(ptr);
#else
  std::free(ptr);
#endif  // defined(_MSC_VER) || defined(__MINGW32__)
}


/*!
 * @brief std::unique_ptr で利用するアラインされたメモリ用のカスタムデリータ
 */
struct AlignedDeleter
{
  /*!
   * @brief デリート処理を行うオペレータ
   * @param [in,out] p  アラインメントされたメモリ領域へのポインタ
   */
  void
  operator()(void* p) const noexcept
  {
    alignedFree(p);
  }
};


/*!
 * @brief プラットフォームIDを取得
 * @param [in] nPlatformEntry  取得するプラットフォームID数の上限
 * @return  プラットフォームIDを格納した std::vector
 */
static inline std::vector<cl_platform_id>
getPlatformIds(cl_uint nPlatformEntry = kNDefaultPlatformEntry)
{
  std::vector<cl_platform_id> platformIds(nPlatformEntry);
  cl_uint nPlatform;
  if (clGetPlatformIDs(nPlatformEntry, platformIds.data(), &nPlatform) != CL_SUCCESS) {
    std::cerr << "clGetPlatformIDs() failed" << std::endl;
    std::exit(EXIT_FAILURE);
  }
  platformIds.resize(nPlatform);
  return platformIds;
}


/*!
 * @brief デバイスIDを取得
 * @param [in] platformId    デバイスIDの取得元のプラットフォームのID
 * @param [in] nDeviceEntry  取得するデバイスID数の上限
 * @param [in] deviceType    取得対象とするデバイス
 * @return デバイスIDを格納した std::vector
 */
static inline std::vector<cl_device_id>
getDeviceIds(const cl_platform_id& platformId, cl_uint nDeviceEntry = kNDefaultDeviceEntry, cl_int deviceType = CL_DEVICE_TYPE_DEFAULT)
{
  std::vector<cl_device_id> deviceIds(nDeviceEntry);
  cl_uint nDevice;
  if (clGetDeviceIDs(platformId, deviceType, nDeviceEntry, deviceIds.data(), &nDevice) != CL_SUCCESS) {
    std::cerr << "clGetDeviceIDs() failed" << std::endl;
    std::exit(EXIT_FAILURE);
  }
  deviceIds.resize(nDevice);
  return deviceIds;
}



/*!
 * @brief カーネル関数へ引数をまとめてセットする関数の実態
 * @param [in] kernel  OpenCLカーネルオブジェクト
 * @param [in] idx     セットする引数のインデックス
 * @param [in] first   セットする引数.可変パラメータから1つだけ取り出したもの
 * @param [in] rest    残りの引数
 * @return OpenCLのエラーコード.エラーが出た時点でエラーコードを返却する.
 */
template<typename First, typename... Rest>
static inline cl_uint
setKernelArgsImpl(const cl_kernel& kernel, int idx, const First& first, const Rest&... rest) noexcept
{
  cl_uint errCode = clSetKernelArg(kernel, idx, sizeof(first), &first);
  return errCode == CL_SUCCESS ? setKernelArgsImpl(kernel, idx + 1, rest...) : errCode;
}


/*!
 * @brief カーネル関数へ最後の引数をセットする
 * @param [in] kernel  OpenCLカーネルオブジェクト
 * @param [in] idx     引数のインデックス
 * @param [in] last    最後の引数
 * @return  OpenCLのエラーコード
 */
template<typename Last>
static inline cl_uint
setKernelArgsImpl(const cl_kernel& kernel, int idx, const Last& last) noexcept
{
  return clSetKernelArg(kernel, idx, sizeof(last), &last);
}


/*!
 * @brief カーネル関数へ引数をまとめてセットする
 * @param [in] kernel  OpenCLカーネルオブジェクト
 * @param [in] args    セットする引数群
 * @return OpenCLのエラーコード.エラーが出た時点でエラーコードを返却する.
 */
template<typename... Args>
static inline cl_uint
setKernelArgs(const cl_kernel& kernel, const Args&... args) noexcept
{
  return setKernelArgsImpl(kernel, 0, args...);
}




/*!
 * @brief このプログラムのエントリポイント
 * @return 終了ステータス
 */
int
main(int argc, char* argv[])
{
  static constexpr int ALIGN = 4096;
  static constexpr std::size_t N = 65536;

  if (argc < 2) {
    std::cerr << "Please specify only one or more source file" << std::endl;
    return EXIT_FAILURE;
  }

  // ホストのメモリを確保
  std::unique_ptr<float[], AlignedDeleter> hostX(alignedMalloc<float*>(N * sizeof(float), ALIGN));
  std::unique_ptr<float[], AlignedDeleter> hostY(alignedMalloc<float*>(N * sizeof(float), ALIGN));
  std::unique_ptr<float[], AlignedDeleter> hostZ(alignedMalloc<float*>(N * sizeof(float), ALIGN));

  // 初期化
  std::mt19937 mt((std::random_device())());
  for (std::size_t i = 0; i < N; i++) {
    hostX[i] = static_cast<float>(mt());
    hostY[i] = static_cast<float>(mt());
  }
  std::fill_n(hostZ.get(), N, 0.0f);

  // プラットフォームを取得
  std::vector<cl_platform_id> platformIds = getPlatformIds(1);

  // デバイスを取得
  std::vector<cl_device_id> deviceIds = getDeviceIds(platformIds[0], 1, CL_DEVICE_TYPE_DEFAULT);

  // コンテキストを生成
  cl_int errCode;
  std::unique_ptr<std::remove_pointer<cl_context>::type, decltype(&clReleaseContext)> context(
      clCreateContext(nullptr, 1, &deviceIds[0], nullptr, nullptr, &errCode), clReleaseContext);

  // コマンドキューを生成
  std::unique_ptr<std::remove_pointer<cl_command_queue>::type, decltype(&clReleaseCommandQueue)> cmdQueue(
      clCreateCommandQueue(context.get(), deviceIds[0], 0, &errCode), clReleaseCommandQueue);

  // デバイスが用いるメモリオブジェクトの生成
  std::unique_ptr<std::remove_pointer<cl_mem>::type, decltype(&clReleaseMemObject)> deviceX(
      clCreateBuffer(context.get(), CL_MEM_READ_WRITE, N * sizeof(float), nullptr, &errCode), clReleaseMemObject);
  std::unique_ptr<std::remove_pointer<cl_mem>::type, decltype(&clReleaseMemObject)> deviceY(
      clCreateBuffer(context.get(), CL_MEM_READ_WRITE, N * sizeof(float), nullptr, &errCode), clReleaseMemObject);
  std::unique_ptr<std::remove_pointer<cl_mem>::type, decltype(&clReleaseMemObject)> deviceZ(
      clCreateBuffer(context.get(), CL_MEM_READ_WRITE, N * sizeof(float), nullptr, &errCode), clReleaseMemObject);

  // ホストのメモリをデバイスのメモリに転送
  errCode = clEnqueueWriteBuffer(cmdQueue.get(), deviceX.get(), CL_TRUE, 0, N * sizeof(float), hostX.get(), 0, nullptr, nullptr);
  errCode = clEnqueueWriteBuffer(cmdQueue.get(), deviceY.get(), CL_TRUE, 0, N * sizeof(float), hostY.get(), 0, nullptr, nullptr);
  errCode = clEnqueueWriteBuffer(cmdQueue.get(), deviceZ.get(), CL_TRUE, 0, N * sizeof(float), hostZ.get(), 0, nullptr, nullptr);

  // コンパイル後のカーネルのバイナリを読み込み
  std::ifstream ifs(argv[1], std::ios::binary);
  if (!ifs.is_open()) {
    std::cerr << "Failed to kernel binary: " << argv[1] << std::endl;
    std::exit(EXIT_FAILURE);
  }
  std::string kernelBin((std::istreambuf_iterator<char>(ifs)), std::istreambuf_iterator<char>());

  // プログラムオブジェクトの生成
  const unsigned char* kbin = reinterpret_cast<const unsigned char*>(kernelBin.c_str());
  std::size_t kbinSize = kernelBin.size();
  cl_int binStatus;
  std::unique_ptr<std::remove_pointer<cl_program>::type, decltype(&clReleaseProgram)> program(
      clCreateProgramWithBinary(context.get(), 1, &deviceIds[0], &kbinSize, &kbin, &binStatus, &errCode), clReleaseProgram);

  // カーネルソースコードのコンパイル (必要な環境もあるらしい?)
  // errCode = clBuildProgram(program.get(), 1, &deviceIds[0], nullptr, nullptr, nullptr);

  // カーネルオブジェクトの生成
  std::unique_ptr<std::remove_pointer<cl_kernel>::type, decltype(&clReleaseKernel)> kernel(
      clCreateKernel(program.get(), "vecAdd", &errCode), clReleaseKernel);

  // カーネル関数に引数を渡す
  errCode = setKernelArgs(kernel.get(), deviceZ.get(), deviceX.get(), deviceY.get(), static_cast<int>(N));

  // カーネルプログラムの実行
  errCode = clEnqueueTask(cmdQueue.get(), kernel.get(), 0, nullptr, nullptr);

  // 終了待機等
  errCode = clFlush(cmdQueue.get());
  errCode = clFinish(cmdQueue.get());

  // 実行結果をデバイスからホストへコピー
  errCode = clEnqueueReadBuffer(cmdQueue.get(), deviceZ.get(), CL_TRUE, 0, N * sizeof(float), hostZ.get(), 0, nullptr, nullptr);

  // 計算結果の確認
  for (std::size_t i = 0; i < N; i++) {
    if (std::abs(hostX[i] + hostY[i] - hostZ[i]) > 1.0e-5) {
      std::cerr << "Result verification failed at element " << i << "!" << std::endl;
      return EXIT_FAILURE;
    }
  }
  std::cout << "Test PASSED" << std::endl;

  return EXIT_SUCCESS;
}

実行は以下のようにする.

$ ./main kernel.bin

これで問題なく実行できればOKである.

雑感

個別のコンパイルプログラムを準備し,事前にオフラインコンパイルを行い,カーネルバイナリを用意するのは正直なところどうなのかと感じた. プログラムの初回起動時に,オンラインコンパイルの結果を保存し,2回目以降にその結果を再利用する形にするのがよさそうだ.

個別のコンパイルプログラムは,OpenCLのプログラムのコンパイルエラーのチェック等に用いるとよいだろう.

もう少しまともな形として作るならば,以下のような形として作るとよさそうだ.

koturn/oclc

まとめ

OpenCLのオフラインコンパイルには,オンラインコンパイルの結果を保存する手法がある.

Osaka.vim #7に行ってきた

Vim

はじめに

参加感想記事なので,今回は「ですます調」で書きます.

2ヶ月遅れの参加記事執筆となりますが,7月2日にOsaka.vimに行ってきました. 7回目のOsaka.vimだったわけですが,6回目の参加でした. ちなみに,Osaka.vimの参加記事を書くのは今回が初めてです(怠惰).

合流

大阪駅に11時頃に到着するとツイートしたら,ryunixさんに「お昼ご飯一緒にどうですか」とお誘いを受けたので,ryunixさん,thincaさん,つじけんさんと合流. Osaka.vim #1 のとき,thincaさんが参加されてたので,何とか合流できました. 会ったことが無い人が多いと,合流に少し苦労しますね.

お昼ご飯

適当に駅周辺をブラつきながら,ご飯食べるところを探してました. 駅前にいたミュージシャン(?)の演奏がうまかった!

で,しばらく歩いてお好み焼き屋に. 店員が怖かったけど,美味しかったです.

本編

今回は大学のときに所属していた研究室の後輩を誘っていたので,隣に座ってVimの知識を叩き込んでいました. 僕は最初からプラグインを導入するより,基本設定である程度やっていた方がよいと考えている派なので,最初は基本的な設定を教えました

とりあえず,カーソルキーを無効化する設定を書かせました.

noremap <Up> <Nop>
noremap <Down> <Nop>
noremap <Right> <Nop>
noremap <Left> <Nop>

最近,僕はVimに行番号は必要ない派になったので,「 set number は不要やで」と言ったのですが,それは受け入れられなかったようですw あとは,ステータスラインを表示する設定とか.

set statusline=%F%m%r%h%w\ [Format=%{&ff}]\ [TYPE=%Y]\ [ASCII=\%03.3b]\ [%p%%]\ {LEN=%L]
set laststatus=2

で,Vim7.3.xxxみたいな化石バージョンのVimを使っていたため,「Vim7.3は古い」という話になり,homebrewで新しいVimを導入させました. それに伴い,プラグインマネージャを導入し,2, 3程度のプラグインを入れるか,という話になり,結局プラグインを導入させる設定を書きました. 日本におけるプラグインマネージャといえば,dein.vimを導入さ,以下のような雛形を渡しました.

if has('vim_starting')
  let s:deindir = expand('~/.cache/dein')
  let s:deinlocal = s:deindir . '/repos/github.com/Shougo/dein.vim'
  let &runtimepath = s:deinlocal . ',' . &runtimepath
endif

if dein#load_state(s:deindir)
  call dein#begin(s:deindir)
  call dein#add('Shougo/dein.vim')
  call dein#add('itchyny/lightline.vim')

  " ここに
  " call dein#add(...)
  " でプラグインを追加する


  call dein#end()
  call dein#save_state()
endif

filetype plugin indent on
syntax enable

(しれっと雛形にサンプルとしてlightline.vimを入れておいたので,後輩のステータスラインがカッコよくなりました)

基本的にはこれでキャッシュ読み込みありでdein.vimを利用できるようになります. TOMLでプラグイン管理も可能なのですが,それは一旦置いておきました.

他には日本語のヘルプを導入させたり.... そんな感じで後輩のVim設定を増やしました.

僕自身の進捗は微妙でしたw やってたことは,job機能をいじるってことで,stdin/stdoutをゴニョゴニョとやってました. ある程度jobの使い方がわかってきたので,これからのプラグイン作成に活かしたいですね.

最後に

やっぱり,こういうVim会に参加すると,Vimに対する意欲が高まるのでいいですね.

Dein.vim で NeoBundleList と同様のことをする

Vim

はじめに

発端はこの会話.

この質問を受け,dein.vimにおいて, :NeoBundleList に相当する関数を探したが,直接的に該当するものはなかった. しかし,dein.vimの関数を組み合わせることで実現できそうだったので,実装を試みた.

実装方針

:NeoBundleList では,管理している全プラグインがそれぞれ

  1. インストールされているかどうか
  2. sourceされている(runtimepath に追加されている)かどうか

echomsg で表示する. これを目指す.

dein#get() を引数を与えずに呼び出すことで,管理プラグインの辞書を得ることができる. この辞書から必要な情報を抜き出すとよい.

インストールされているかどうか

dein.vimが想定するインストール先のパスを確認するとよい. イメージとしては以下の通り.

for pair in items(dein#get())
  echomsg pair[0] 'is' (isdirectory(pair[1].path) ? 'installed' : 'not installed')
endfor

sourceされているかどうか

それぞれのプラグインの辞書に sourced というキーがあるので,それを利用する. イメージとしては以下の通り.

for pair in items(dein#get())
  echomsg pair[0] 'is' (pair[1].sourced ? 'sourced' : 'not sourced')
endfor

実装

以上のことをまとめ,以下のコードにするとよい. ついでに, DeinList というコマンドを定義しておく.

function! s:dein_list() abort
  echomsg '[dein] #: not sourced, X: not installed'
  for pair in items(dein#get())
    echomsg (!isdirectory(pair[1].path) ? 'X'
          \ : pair[1].sourced ? ' '
          \ : '#') pair[0]
  endfor
endfunction
command! DeinList call s:dein_list()

ちなみに,空引数の dein#get()dein#_plugins の浅いコピーを返却するだけなので,以下のようにした方が無駄がない.

function! s:dein_list() abort
  echomsg '[dein] #: not sourced, X: not installed'
  for pair in items(dein#_plugins)
    echomsg (!isdirectory(pair[1].path) ? 'X'
          \ : pair[1].sourced ? ' '
          \ : '#') pair[0]
  endfor
endfunction
command! DeinList call s:dein_list()

ただ,privateを意図したような変数を直接いじるのは,やや気乗りしないものがある.

まとめ

dein.vimの関数を組み合わせることで,:NeoBundleList と同等のものを実現することができる.

MSYS2でOpenCLを使う

C++ OpenCL

はじめに

MSYS2でOpenCLのプログラムをコンパイル&実行したかった. pacmanOpenCL関連のヘッダを導入することはできるが,OpenCLのインポートライブラリは導入することはできない. ここでは,MSYS2でOpenCLの環境を導入する一連の手順を紹介する.

OpenCLのヘッダを導入する

pacmanを用いると,OpenCLのヘッダを導入できる. x64ならば,

$ pacman -S mingw-w64-x86_64-opencl-headers

x86ならば,

$ pacman -S mingw-w64-i686-opencl-headers

とするとよい.

OpenCLのインポートライブラリを作成する

上記の手順では,ヘッダファイルしか導入できない. インポートライブラリは自分で作成する必要がある. この手順は,ここを参考にした.

  • x64環境
$ mkdir lib64 && gendef - /c/Windows/system32/OpenCL.dll > lib64/OpenCL.def
$ dlltool -l lib64/libOpenCL.a -d lib64/OpenCL.def -A -k
$ mkdir lib32 && gendef - /c/Windows/SysWOW64/OpenCL.dll > lib32/OpenCL.def
$ dlltool -l lib32/libOpenCL.a -d lib32/OpenCL.def -A -k

あとは必要に応じて, /usr/local/lib 下に置くなどするとよい.

コンパイルの確認

適当に以下のようなソースコードを準備する. (面倒なので, main() 関数でのエラー処理は一切行っていない.)

#include <cstdlib>
#include <fstream>
#include <iostream>
#include <memory>
#include <random>
#include <string>

#ifdef __APPLE__
#  include <OpenCL/opencl.h>
#else
#  include <CL/cl.h>
#endif


/*!
 * @brief アラインメントされたメモリを動的確保する関数
 * @param [in] size       確保するメモリサイズ (単位はbyte)
 * @param [in] alignment  アラインメント (2のべき乗を指定すること)
 * @return  アラインメントし,動的確保されたメモリ領域へのポインタ
 */
template<typename T = void*, typename std::enable_if<std::is_pointer<T>::value, std::nullptr_t>::type = nullptr>
static inline T
alignedMalloc(std::size_t size, std::size_t alignment) noexcept
{
#if defined(_MSC_VER) || defined(__MINGW32__)
  return reinterpret_cast<T>(_aligned_malloc(size, alignment));
#else
  void* p;
  return reinterpret_cast<T>(posix_memalign(&p, alignment, size) == 0 ? p : nullptr);
#endif  // defined(_MSC_VER) || defined(__MINGW32__)
}


/*!
 * @brief アラインメントされたメモリを解放する関数
 * @param [in] ptr  解放対象のメモリの先頭番地を指すポインタ
 */
static inline void
alignedFree(void* ptr) noexcept
{
#if defined(_MSC_VER) || defined(__MINGW32__)
  _aligned_free(ptr);
#else
  std::free(ptr);
#endif  // defined(_MSC_VER) || defined(__MINGW32__)
}


/*!
 * @brief std::unique_ptr で利用するアラインされたメモリ用のカスタムデリータ
 */
struct AlignedDeleter
{
  /*!
   * @brief デリート処理を行うオペレータ
   * @param [in,out] p  アラインメントされたメモリ領域へのポインタ
   */
  void
  operator()(void* p) const noexcept
  {
    alignedFree(p);
  }
};


/*!
 * @brief OpenCLのデバイス情報を表示する
 * @param [in] device  OpenCLのデバイスID
 */
static inline void
showDeviceInfo(cl_device_id device) noexcept
{
  char info[2048];
  std::size_t size;

  int err = clGetDeviceInfo(device, CL_DEVICE_NAME, sizeof(info), info, &size);
  std::cout << "  Device: " << info << std::endl;
  err = clGetDeviceInfo(device, CL_DEVICE_VERSION, sizeof(info), info, &size);
  std::cout << "  CL_DEVICE_VERSION: " << info << std::endl;
}


/*!
 * @brief OpenCLのプラットフォーム情報を表示する
 * @param [in] platformId  OpenCLのプラットフォーム情報
 */
static inline void
showPlatformInfo(cl_platform_id platformId) noexcept
{
  static constexpr int MAX_DEVICE_IDS = 8;
  char info[2048];

  std::cout << "========== Platform Information ==========" << std::endl;
  std::size_t size;
  int err = clGetPlatformInfo(platformId, CL_PLATFORM_NAME, sizeof(info), info, &size);
  std::cout << "Platform: " << info << std::endl;
  err = clGetPlatformInfo(platformId, CL_PLATFORM_VERSION, sizeof(info), info, &size);
  std::cout << "CL_PLATFORM_VERSION: " << info << std::endl;

  cl_device_id devices[MAX_DEVICE_IDS]; cl_uint nDevice;
  err = clGetDeviceIDs(platformId, CL_DEVICE_TYPE_ALL, MAX_DEVICE_IDS, devices, &nDevice);
  if (err != CL_SUCCESS) {
    std::cerr << "Error (clGetDeviceIDs): " << err << std::endl;
    std::exit(EXIT_FAILURE);
  }
  for (cl_uint j = 0; j < nDevice; j++) {
    showDeviceInfo(devices[j]);
  }
  std::cout << "==========================================" << std::endl;
}


/*!
 * @brief カーネル関数へ引数をまとめてセットする関数の実態
 * @param [in] kernel  OpenCLカーネルオブジェクト
 * @param [in] idx     セットする引数のインデックス
 * @param [in] first   セットする引数.可変パラメータから1つだけ取り出したもの
 * @param [in] rest    残りの引数
 * @return OpenCLのエラーコード.エラーが出た時点でエラーコードを返却する.
 */
template<typename First, typename... Rest>
static inline cl_uint
setKernelArgsImpl(const cl_kernel& kernel, int idx, const First& first, const Rest&... rest) noexcept
{
  cl_uint errCode = clSetKernelArg(kernel, idx, sizeof(first), &first);
  return errCode == CL_SUCCESS ? setKernelArgsImpl(kernel, idx + 1, rest...) : errCode;
}


/*!
 * @brief カーネル関数へ最後の引数をセットする
 * @param [in] kernel  OpenCLカーネルオブジェクト
 * @param [in] idx     引数のインデックス
 * @param [in] last    最後の引数
 * @return  OpenCLのエラーコード
 */
template<typename Last>
static inline cl_uint
setKernelArgsImpl(const cl_kernel& kernel, int idx, const Last& last) noexcept
{
  return clSetKernelArg(kernel, idx, sizeof(last), &last);
}


/*!
 * @brief カーネル関数へ引数をまとめてセットする
 * @param [in] kernel  OpenCLカーネルオブジェクト
 * @param [in] args    セットする引数群
 * @return OpenCLのエラーコード.エラーが出た時点でエラーコードを返却する.
 */
template<typename... Args>
static inline cl_uint
setKernelArgs(const cl_kernel& kernel, const Args&... args) noexcept
{
  return setKernelArgsImpl(kernel, 0, args...);
}




/*!
 * @brief このプログラムのエントリポイント
 * @return 終了ステータス
 */
int
main()
{
  static constexpr int ALIGN = 4096;
  static constexpr std::size_t N = 65536;
  static const char KERNEL_FILENAME[] = "kernel.cl";

  // ホストのメモリを確保
  std::unique_ptr<float[], AlignedDeleter> hostX(alignedMalloc<float*>(N * sizeof(float), ALIGN));
  std::unique_ptr<float[], AlignedDeleter> hostY(alignedMalloc<float*>(N * sizeof(float), ALIGN));
  std::unique_ptr<float[], AlignedDeleter> hostZ(alignedMalloc<float*>(N * sizeof(float), ALIGN));

  // 初期化
  std::mt19937 mt((std::random_device())());
  for (std::size_t i = 0; i < N; i++) {
    hostX[i] = static_cast<float>(mt());
    hostY[i] = static_cast<float>(mt());
  }
  std::fill_n(hostZ.get(), N, 0.0f);

  // プラットフォームを取得
  cl_platform_id platformId;
  cl_uint nPlatform;
  cl_int errCode = clGetPlatformIDs(1, &platformId, &nPlatform);
  showPlatformInfo(platformId);

  // デバイス情報を取得
  cl_device_id deviceId;
  cl_uint nDevice;
  errCode = clGetDeviceIDs(platformId, CL_DEVICE_TYPE_DEFAULT, 1, &deviceId, &nDevice);

  // コンテキストを生成
  std::unique_ptr<std::remove_pointer<cl_context>::type, decltype(&clReleaseContext)> context(
      clCreateContext(nullptr, 1, &deviceId, nullptr, nullptr, &errCode), clReleaseContext);

  // コマンドキューを生成
  std::unique_ptr<std::remove_pointer<cl_command_queue>::type, decltype(&clReleaseCommandQueue)> cmdQueue(
      clCreateCommandQueue(context.get(), deviceId, 0, &errCode), clReleaseCommandQueue);

  // デバイスが用いるメモリオブジェクトの生成
  std::unique_ptr<std::remove_pointer<cl_mem>::type, decltype(&clReleaseMemObject)> deviceX(
      clCreateBuffer(context.get(), CL_MEM_READ_WRITE, N * sizeof(float), nullptr, &errCode), clReleaseMemObject);
  std::unique_ptr<std::remove_pointer<cl_mem>::type, decltype(&clReleaseMemObject)> deviceY(
      clCreateBuffer(context.get(), CL_MEM_READ_WRITE, N * sizeof(float), nullptr, &errCode), clReleaseMemObject);
  std::unique_ptr<std::remove_pointer<cl_mem>::type, decltype(&clReleaseMemObject)> deviceZ(
      clCreateBuffer(context.get(), CL_MEM_READ_WRITE, N * sizeof(float), nullptr, &errCode), clReleaseMemObject);

  // ホストのメモリをデバイスのメモリに転送
  errCode = clEnqueueWriteBuffer(cmdQueue.get(), deviceX.get(), CL_TRUE, 0, N * sizeof(float), hostX.get(), 0, nullptr, nullptr);
  errCode = clEnqueueWriteBuffer(cmdQueue.get(), deviceY.get(), CL_TRUE, 0, N * sizeof(float), hostY.get(), 0, nullptr, nullptr);
  errCode = clEnqueueWriteBuffer(cmdQueue.get(), deviceZ.get(), CL_TRUE, 0, N * sizeof(float), hostZ.get(), 0, nullptr, nullptr);

  // カーネルのソースコードを文字列として取得
  std::ifstream ifs(KERNEL_FILENAME);
  if (!ifs.is_open()) {
    std::cerr << "Failed to read kernel program: " << KERNEL_FILENAME << std::endl;
    return EXIT_FAILURE;
  }
  std::string kernelSource((std::istreambuf_iterator<char>(ifs)), std::istreambuf_iterator<char>());

  // プログラムオブジェクトの生成
  const char* ksrc = kernelSource.c_str();
  std::string::size_type srcSize = kernelSource.length();
  std::unique_ptr<std::remove_pointer<cl_program>::type, decltype(&clReleaseProgram)> program(
      clCreateProgramWithSource(context.get(), 1, &ksrc, &srcSize, &errCode), clReleaseProgram);

  // カーネルソースコードのコンパイル
  errCode = clBuildProgram(program.get(), 1, &deviceId, nullptr, nullptr, nullptr);

  // カーネルオブジェクトの生成
  std::unique_ptr<std::remove_pointer<cl_kernel>::type, decltype(&clReleaseKernel)> kernel(
      clCreateKernel(program.get(), "vecAdd", &errCode), clReleaseKernel);

  // カーネル関数に引数を渡す
  errCode = setKernelArgs(kernel.get(), deviceZ.get(), deviceX.get(), deviceY.get(), static_cast<int>(N));

  // カーネルプログラムの実行
  errCode = clEnqueueTask(cmdQueue.get(), kernel.get(), 0, nullptr, nullptr);

  // 終了待機等
  errCode = clFlush(cmdQueue.get());
  errCode = clFinish(cmdQueue.get());

  // 実行結果をデバイスからホストへコピー
  errCode = clEnqueueReadBuffer(cmdQueue.get(), deviceZ.get(), CL_TRUE, 0, N * sizeof(float), hostZ.get(), 0, nullptr, nullptr);

  // 計算結果の確認
  for (std::size_t i = 0; i < N; i++) {
    if (std::abs(hostX[i] + hostY[i] - hostZ[i]) > 1.0e-5) {
      std::cerr << "Result verification failed at element " << i << "!" << std::endl;
      return EXIT_FAILURE;
    }
  }
  std::cout << "Test PASSED" << std::endl;

  return EXIT_SUCCESS;
}

そして,カーネルソースコードとして以下のものを準備する. ファイル名は kernel.cl とする.

__kernel void
vecAdd(__global float* z, __global float* x, __global float* y, int n)
{
  const int para = 4;
  const int end = (n / para) * para;

  for (int i = 0; i < end; i += para) {
    float4 vtmp = vload4(0, x + i) + vload4(0, y + i);
    vstore4(vtmp, 0, z + i);
  }

  for (int i = end; i < n; i++) {
    z[i] = x[i] + y[i];
  }
}

コンパイルと実行は以下の通り. これで問題なくコンパイル&実行できればOKである.

$ g++ -std=gnu++11 main.cpp -Llib64 -lOpenCL -o main

余談

カーネル関数へ引数を渡す関数: clSetKernelArg() は,引数のインデックスと引数へのポインタを渡す形となっており,非常にダサい. また,カーネル関数の引数の数だけ呼び出さなければならないのも大変だ.

int n = static_cast<int>(N);
clSetKernelArg(kernel.get(), 0, sizeof(deviceZ.get()), &deviceZ.get());
clSetKernelArg(kernel.get(), 1, sizeof(deviceX.get()), &deviceX.get());
clSetKernelArg(kernel.get(), 2, sizeof(deviceY.get()), &deviceY.get());
clSetKernelArg(kernel.get(), 3, sizeof(n), &n);

うまく可変引数テンプレートを用いれば,このダサさを解消できると考え,以下のような関数の実装を行った.

template<typename First, typename... Rest>
static inline cl_uint
setKernelArgsImpl(const cl_kernel& kernel, int idx, const First& first, const Rest&... rest) noexcept
{
  cl_uint errCode = clSetKernelArg(kernel, idx, sizeof(first), &first);
  return errCode == CL_SUCCESS ? setKernelArgsImpl(kernel, idx + 1, rest...) : errCode;
}


template<typename Last>
static inline cl_uint
setKernelArgsImpl(const cl_kernel& kernel, int idx, const Last& last) noexcept
{
  return clSetKernelArg(kernel, idx, sizeof(last), &last);
}


template<typename... Args>
static inline cl_uint
setKernelArgs(const cl_kernel& kernel, const Args&... args) noexcept
{
  return setKernelArgsImpl(kernel, 0, args...);
}

実装としては,ユーザが呼び出しを行う setKernelArgs() と実際の処理を担当する setKernelArgsImpl() の2つに分けている. setKernelArgsImpl() は,可変テンプレートのアンパックとインデックスのインクリメントを行う. 引数のサイズは sizeof 演算子で取得するので,引数の型はちゃんと意識する必要がある.

呼び出し側は以下のようになる. 呼び出し側の見た目としては,

  • インデックスが必要ない
  • カーネル関数の引数をポインタとして渡す必要がない
  • 引数のサイズを渡す必要がない

ため,非常にスッキリしていると思う.

setKernelArgs(kernel.get(), deviceZ.get(), deviceX.get(), deviceY.get(), static_cast<int>(N));

C++11様々である.

まとめ

MSYS2でOpenCLを使用するには,pacmanでOepnCL関連のヘッダのインストールを行い,WindowsOpenCLのDLLからインポートライブラリを作成する必要がある.

参考文献

std::vectorなどでアラインされた領域を用いる

C++ Qt

はじめに

前回の記事でCPU SIMD命令(SSE/AVX/NEON)を紹介した. ただ,前回の記事では静的配列や単純な動的確保のパターンでしかアラインされた領域を用いる手法しか紹介しなかった.

今回は,C++でよく用いられるであろう std::vector でアラインされた領域を用いる手法を紹介する. これにより, std::vector などのSTLのデータにもSIMD処理を適用することが可能となる. また, std::unique_ptrstd::shared_ptr でアラインされた領域を解放する手法も紹介する.

STLのカスタムアロケータ

std::vector などのSTLメモリ確保は基本的に ::operator new を用いて行われる. これのメモリ確保は std::allocator というクラスが担当し, STL内部で用いられている.

例えば, std::vector のテンプレートの第二引数を省略した場合, std::allocator が用いられる. この std::vector のテンプレート第二引数に自作アロケータを指定すると, std::vector でもアラインされたメモリを使用できる.

さて, std::vector のアロケータの自作であるが,以下を参考にした. std::allocator の実装が簡略化されて紹介されている.

カスタムアロケータの自作は std::allocator を継承する形で実装する. オーバーライドすべきメンバ関数は,

  • allocate()
  • deallocate()

の2つである. 実装例としては以下のようになる.

// $ g++ -std=gnu++11 main.cpp -o main
#include <cstdlib>
#include <iostream>
#include <type_traits>
#include <vector>


/*!
 * @brief アラインメントされたメモリを動的確保する関数
 * @param [in] size       確保するメモリサイズ (単位はbyte)
 * @param [in] alignment  アラインメント (2のべき乗を指定すること)
 * @return  アラインメントし,動的確保されたメモリ領域へのポインタ
 */
template<typename T = void*, typename std::enable_if<std::is_pointer<T>::value, std::nullptr_t>::type = nullptr>
static inline T
alignedMalloc(std::size_t size, std::size_t alignment) noexcept
{
#if defined(_MSC_VER) || defined(__MINGW32__)
  return reinterpret_cast<T>(_aligned_malloc(size, alignment));
#else
  void* p;
  return reinterpret_cast<T>(posix_memalign(&p, alignment, size) == 0 ? p : nullptr);
#endif  // defined(_MSC_VER) || defined(__MINGW32__)
}


/*!
 * @brief アラインメントされたメモリを解放する関数
 * @param [in] ptr  解放対象のメモリの先頭番地を指すポインタ
 */
static inline void
alignedFree(void* ptr) noexcept
{
#if defined(_MSC_VER) || defined(__MINGW32__)
  _aligned_free(ptr);
#else
  std::free(ptr);
#endif  // defined(_MSC_VER) || defined(__MINGW32__)
}


/*!
 * @brief アラインメントされたメモリを確保するSTLのカスタムアロケータ
 */
template<typename T, std::size_t N>
class AligndSTLAllocator :
  public std::allocator<T>
{
public:
  using size_type = typename std::allocator<T>::size_type;
  using pointer = typename std::allocator<T>::pointer;
  using const_pointer = typename std::allocator<T>::const_pointer;

  /*!
   * @brief STLのメモリ領域を動的確保する
   *
   * @param [in] n              確保する要素数
   * @param [in] const_pointer  使用しない引数
   * @return  アラインされたメモリ領域へのポインタ
   */
  pointer
  allocate(size_type n, const_pointer = nullptr) const
  {
    if (n > this->max_size()) {
      throw std::bad_alloc("Cannot allocate aligned memory.");
    }
    return alignedMalloc<pointer>(n * sizeof(T), N);
  }

  /*!
   * @brief STLのメモリ領域を解放する
   * @param [in,out] p  動的確保したメモリ領域
   * @param [in]     n  要素数 (使用しない引数)
   */
  void
  deallocate(pointer p, size_type) const noexcept
  {
    alignedFree(p);
  }
};


/*!
 * @brief このプログラムのエントリポイント
 * @return 終了ステータス
 */
int
main()
{
  static constexpr int ALIGN = 32;

  std::vector<int, AligndSTLAllocator<int, ALIGN> > vec(32);
  if ((reinterpret_cast<std::ptrdiff_t>(vec.data())) % ALIGN == 0) {
    std::cout << "vector memory is " << ALIGN << " byte aligned." << std::endl;
  } else {
    std::cout << "vector memory is not " << ALIGN << " byte aligned." << std::endl;
  }

  return EXIT_SUCCESS;
}

std::unique_ptrstd::shared_ptr のカスタムデリータ

std::unique_ptrstd::shared_ptr はポインタを管理するクラスであり,自動的に始末処理を行ってくれる. 始末処理はユーザが指定することが可能である. 指定の方法として,

  • 関数オブジェクト
  • 関数ポインタ
  • ラムダ

を用いるものの3つがある. これらの実装例は以下のようになる.

// $ g++ -std=gnu++11 main.cpp -o main
#include <cstdlib>
#include <iostream>
#include <memory>
#include <type_traits>


/*!
 * @brief アラインメントされたメモリを動的確保する関数
 * @param [in] size       確保するメモリサイズ (単位はbyte)
 * @param [in] alignment  アラインメント (2のべき乗を指定すること)
 * @return  アラインメントし,動的確保されたメモリ領域へのポインタ
 */
template<typename T = void*, typename std::enable_if<std::is_pointer<T>::value, std::nullptr_t>::type = nullptr>
static inline T
alignedMalloc(std::size_t size, std::size_t alignment) noexcept
{
#if defined(_MSC_VER) || defined(__MINGW32__)
  return reinterpret_cast<T>(_aligned_malloc(size, alignment));
#else
  void* p;
  return reinterpret_cast<T>(posix_memalign(&p, alignment, size) == 0 ? p : nullptr);
#endif  // defined(_MSC_VER) || defined(__MINGW32__)
}


/*!
 * @brief アラインメントされたメモリを解放する関数
 * @param [in] ptr  解放対象のメモリの先頭番地を指すポインタ
 */
static inline void
alignedFree(void* ptr) noexcept
{
#if defined(_MSC_VER) || defined(__MINGW32__)
  _aligned_free(ptr);
#else
  std::free(ptr);
#endif  // defined(_MSC_VER) || defined(__MINGW32__)
}


/*!
 * @brief std::unique_ptr で利用するアラインされたメモリ用のカスタムデリータ
 */
struct AlignedDeleter
{
  /*!
   * @brief デリート処理を行うオペレータ
   * @param [in,out] p  アラインメントされたメモリ領域へのポインタ
   */
  void
  operator()(void* p) const noexcept
  {
    alignedFree(p);
  }
};


/*!
 * @brief アラインメントをチェックする
 * @param [in] name       チェック対象ポインタ名
 * @param [in] ptr        チェック対象ポインタ
 * @param [in] alignment  期待するアラインメント
 */
static inline void
checkAlignment(const std::string& name, const void* ptr, std::size_t alignment) noexcept
{
  if ((reinterpret_cast<std::ptrdiff_t>(ptr)) % alignment == 0) {
    std::cout << name << " is " << alignment << " byte aligned." << std::endl;
  } else {
    std::cout << name << " is not " << alignment << " byte aligned." << std::endl;
  }
}


/*!
 * @brief このプログラムのエントリポイント
 * @return 終了ステータス
 */
int
main()
{
  static constexpr int ALIGN = 32;
  static constexpr int N = 10;

  // std::unique_ptr にカスタムデリータを指定する例1: 関数オブジェクト
  std::unique_ptr<unsigned char, AlignedDeleter> up01(alignedMalloc<unsigned char*>(N * sizeof(unsigned char), ALIGN));
  if (up01.get() == nullptr) {
    std::cerr << "Failed to allocate memory." << std::endl;
    return 1;
  }
  checkAlignment("up01", up01.get(), ALIGN);

  // std::unique_ptr にカスタムデリータを指定する例2: 関数ポインタ
  std::unique_ptr<unsigned char, decltype(&alignedFree)> up02(alignedMalloc<unsigned char*>(N * sizeof(unsigned char), ALIGN), alignedFree);
  if (up02.get() == nullptr) {
    std::cerr << "Failed to allocate memory." << std::endl;
    return 1;
  }
  checkAlignment("up02", up02.get(), ALIGN);

  // std::unique_ptr にカスタムデリータを指定する例3: ラムダ
  auto deleter = [](void* ptr){
    alignedFree(ptr);
  };
  std::unique_ptr<unsigned char, decltype(deleter)> up03(alignedMalloc<unsigned char*>(N * sizeof(unsigned char), ALIGN), deleter);
  if (up03.get() == nullptr) {
    std::cerr << "Failed to allocate memory." << std::endl;
    return 1;
  }
  checkAlignment("up03", up03.get(), ALIGN);

  // std::shared_ptr にカスタムデリータを指定する例1: 関数オブジェクト
  std::shared_ptr<unsigned char> sp01(alignedMalloc<unsigned char*>(N * sizeof(unsigned char), ALIGN), AlignedDeleter());
  if (sp01.get() == nullptr) {
    std::cerr << "Failed to allocate memory." << std::endl;
    return 1;
  }
  checkAlignment("sp01", sp01.get(), ALIGN);

  // std::shared_ptr にカスタムデリータを指定する例2: 関数ポインタ
  std::shared_ptr<unsigned char> sp02(alignedMalloc<unsigned char*>(N * sizeof(unsigned char), ALIGN), alignedFree);
  if (sp02.get() == nullptr) {
    std::cerr << "Failed to allocate memory." << std::endl;
    return 1;
  }
  checkAlignment("sp02", sp02.get(), ALIGN);

  // std::shared_ptr にカスタムデリータを指定する例3: ラムダ
  std::shared_ptr<unsigned char> sp03(alignedMalloc<unsigned char*>(N * sizeof(unsigned char), ALIGN), [](void* ptr) {
    alignedFree(ptr);
  });
  if (sp03.get() == nullptr) {
    std::cerr << "Failed to allocate memory." << std::endl;
    return 1;
  }
  checkAlignment("sp03", sp03.get(), ALIGN);

  return 0;
}

おまけ: Qtの QScopedPointerQSharedPointer のカスタムデリータ

Qtには std::unique_ptr に相当する QScopedPointerstd::shared_ptr に相当する QSharedPointer がある. これらのカスタムデリータについてもおまけで紹介する.

#include <QCoreApplication>
#include <QDebug>
#include <QScopedPointer>
#include <QScopedArrayPointer>
#include <QSharedPointer>

#include <string>


/*!
 * @brief アラインメントされたメモリを動的確保する関数
 * @param [in] size       確保するメモリサイズ (単位はbyte)
 * @param [in] alignment  アラインメント (2のべき乗を指定すること)
 * @return  アラインメントし,動的確保されたメモリ領域へのポインタ
 */
template<typename T = void*, typename std::enable_if<std::is_pointer<T>::value, std::nullptr_t>::type = nullptr>
static inline T
alignedMalloc(std::size_t size, std::size_t alignment) noexcept
{
#if defined(_MSC_VER) || defined(__MINGW32__)
  return reinterpret_cast<T>(_aligned_malloc(size, alignment));
#else
  void* p;
  return reinterpret_cast<T>(posix_memalign(&p, alignment, size) == 0 ? p : nullptr);
#endif  // defined(_MSC_VER) || defined(__MINGW32__)
}


/*!
 * @brief アラインメントされたメモリを解放する関数
 * @param [in] ptr  解放対象のメモリの先頭番地を指すポインタ
 */
static inline void
alignedFree(void* ptr) noexcept
{
#if defined(_MSC_VER) || defined(__MINGW32__)
  _aligned_free(ptr);
#else
  std::free(ptr);
#endif  // defined(_MSC_VER) || defined(__MINGW32__)
}


/*!
 * @brief QScopedPointer, QScopedArrayPointer, QSharedPointer の
 *        カスタムデリータを提供する
 */
struct AlignedQDeleter
{
  /*!
   * @brief QScopedPointer, QScopedArrayPointer で利用するデリート関数
   * @param [in,out] ptr  アラインメントされたメモリ領域へのポインタ
   */
  static inline void
  cleanup(void* ptr) noexcept
  {
    alignedFree(ptr);
  }

  /*!
   * @brief QSharedPointer で利用するデリート関数
   * @param [in,out] ptr  アラインメントされたメモリ領域へのポインタ
   */
  inline void
  operator()(void* ptr) const noexcept
  {
    cleanup(ptr);
  }
};


/*!
 * @brief アラインメントをチェックする
 * @param [in] name       チェック対象ポインタ名
 * @param [in] ptr        チェック対象ポインタ
 * @param [in] alignment  期待するアラインメント
 */
static inline void
checkAlignment(const std::string& name, const void* ptr, std::size_t alignment) noexcept
{
  if ((reinterpret_cast<std::ptrdiff_t>(ptr)) % alignment == 0) {
    qDebug() << name.c_str() << "is" << alignment << "byte aligned.";
  } else {
    qDebug() << name.c_str() << "is not" << alignment << "byte aligned.";
  }
}


/*!
 * @brief このプログラムのエントリポイント
 * @return 終了ステータス
 */
int
main(int argc, char *argv[])
{
  static constexpr int ALIGN = 32;
  static constexpr int N = 10;

  QCoreApplication a(argc, argv);

  QScopedPointer<quint8, AlignedQDeleter> scopedPtr(alignedMalloc<quint8*>(N * sizeof(quint8), ALIGN));
  if (scopedPtr.data() == nullptr) {
    qDebug() << "Failed to allocate memory.";
    return 1;
  }
  checkAlignment("scopedPtr", scopedPtr.data(), ALIGN);

  QScopedArrayPointer<quint8, AlignedQDeleter> scopedArrayPtr(alignedMalloc<quint8*>(N * sizeof(quint8), ALIGN));
  if (scopedArrayPtr.data() == nullptr) {
    qDebug() << "Failed to allocate memory.";
    return 1;
  }
  checkAlignment("scopedArrayPtr", scopedArrayPtr.data(), ALIGN);

  QSharedPointer<quint8> sharedPtr(alignedMalloc<quint8*>(N * sizeof(quint8), ALIGN), AlignedQDeleter());
  if (sharedPtr.data() == nullptr) {
    qDebug() << "Failed to allocate memory.";
    return 1;
  }
  checkAlignment("sharedPtr", sharedPtr.data(), ALIGN);

  return a.exec();
}

コンパイルするには,以下のようなプロジェクトファイルを準備する.

QT += core
QT -= gui

CONFIG += c++11

TARGET = CustomDeleterTest
CONFIG += console
CONFIG -= app_bundle

TEMPLATE = app

SOURCES += main.cpp

そして,以下の手順でコンパイルを行うとよい.

$ qmake
$ make

まとめ

std::vector にカスタムアロケータを指定することで,データのアラインメントを指定できる. これにより, std::vector SIMD処理に用いることが可能となる.

また. std::unique_ptrstd::shared_ptr でアラインされた領域を管理する場合,カスタムデリータで解放処理を指定する必要がある.

参考文献

AVXとAVX-512のインタリーブ

C++

はじめに

前回の記事では,Intel系のCPUとARM系のCPUのSIMD命令紹介した. 記事中のサンプルコードで,画像の2倍の拡大を行うコードがあり,その中でインタリーブを行っていた.

SSEであれば,単純にunpack命令を実行するだけでよかった. 簡単なサンプルコードと出力結果は以下の通り.

// $ g++ -std=gnu++14 -march=native main.cpp -o main
#include <iostream>
#ifdef _MSC_VER
#  include <intrin.h>
#else
#  include <x86intrin.h>
#endif  // _MSC_VER


static inline void
showAVX128(const __m128i& v128) noexcept
{
  alignas(alignof(__m128i)) unsigned char z[sizeof(__m128i)] = {0};
  _mm_store_si128(reinterpret_cast<__m128i*>(z), v128);
  for (const auto& b : z) {
    std::cout << static_cast<unsigned int>(b) << ' ';
  }
  std::cout << "\n" << std::endl;
}


int
main()
{
  alignas(alignof(__m128i)) unsigned char z[sizeof(__m128i)] = {0};

  for (std::size_t i = 0; i < sizeof(z) / sizeof(z[0]); i++) {
    z[i] = static_cast<unsigned char>(i);
  }

  __m128i v128 = _mm_load_si128(reinterpret_cast<const __m128i*>(z));

  std::cout << "==================== low ====================" << std::endl;
  showAVX128(_mm_unpacklo_epi8(v128, v128));

  std::cout << "==================== high ====================" << std::endl;
  showAVX128(_mm_unpackhi_epi8(v128, v128));

  return 0;
}
==================== low ====================
0 0 1 1 2 2 3 3 4 4 5 5 6 6 7 7

==================== high ====================
8 8 9 9 10 10 11 11 12 12 13 13 14 14 15 15

しかし,AVX,AVX-512のインタリーブは一筋縄ではいかない. 例えば,以下の2つコード,

// $ g++ -std=gnu++14 -march=native main.cpp -o main
#include <iostream>
#ifdef _MSC_VER
#  include <intrin.h>
#else
#  include <x86intrin.h>
#endif  // _MSC_VER


static inline void
showAVX256(const __m256i& v256) noexcept
{
  alignas(alignof(__m256i)) unsigned char z[sizeof(__m256i)] = {0};
  _mm256_store_si256(reinterpret_cast<__m256i*>(z), v256);
  for (const auto& b : z) {
    std::cout << static_cast<unsigned int>(b) << ' ';
  }
  std::cout << "\n" << std::endl;
}


int
main()
{
  alignas(alignof(__m256i)) unsigned char z[sizeof(__m256i)] = {0};

  for (std::size_t i = 0; i < sizeof(z) / sizeof(z[0]); i++) {
    z[i] = static_cast<unsigned char>(i);
  }

  __m256i v256 = _mm256_load_si256(reinterpret_cast<const __m256i*>(z));

  std::cout << "==================== low ====================" << std::endl;
  showAVX256(_mm256_unpacklo_epi8(v256, v256));

  std::cout << "==================== high ====================" << std::endl;
  showAVX256(_mm256_unpackhi_epi8(v256, v256));

  return 0;
}
// $ g++ -std=gnu++14 -march=native -mavx512vbmi main.cpp -o main
#include <iostream>
#ifdef _MSC_VER
#  include <intrin.h>
#else
#  include <x86intrin.h>
#endif  // _MSC_VER


static inline void
showAVX512(const __m512i& v512) noexcept
{
  alignas(alignof(__m512i)) unsigned char z[sizeof(__m512i)] = {0};
  _mm512_store_si512(reinterpret_cast<__m512i*>(z), v512);
  for (const auto& b : z) {
    std::cout << static_cast<unsigned int>(b) << ' ';
  }
  std::cout << "\n" << std::endl;
}


int
main()
{
  alignas(alignof(__m512i)) unsigned char z[sizeof(__m512i)] = {0};

  for (std::size_t i = 0; i < sizeof(z) / sizeof(z[0]); i++) {
    z[i] = static_cast<unsigned char>(i);
  }

  __m512i v512 = _mm512_load_si512(reinterpret_cast<const __m512i*>(z));

  std::cout << "==================== low ====================" << std::endl;
  showAVX512(_mm512_unpacklo_epi8(v512, v512));

  std::cout << "==================== high ====================" << std::endl;
  showAVX512(_mm512_unpackhi_epi8(v512, v512));

  return 0;
}

の出力はそれぞれ次のようなものを期待する.

==================== low ====================
0 0 1 1 2 2 3 3 4 4 5 5 6 6 7 7 8 8 9 9 10 10 11 11 12 12 13 13 14 14 15 15

==================== high ====================
16 16 17 17 18 18 19 19 20 20 21 21 22 22 23 23 24 24 25 25 26 26 27 27 28 28 29 29 30 30 31 31
==================== low ====================
0 0 1 1 2 2 3 3 4 4 5 5 6 6 7 7 8 8 9 9 10 10 11 11 12 12 13 13 14 14 15 15 16 16 17 17 18 18 19 19 20 20 21 21 22 22 23 23 24 24 25 25 26 26 27 27 28 28 29 29 30 30 31 31

==================== high ====================
32 32 33 33 34 34 35 35 36 36 37 37 38 38 39 39 40 40 41 41 42 42 43 43 44 44 45 45 46 46 47 47 48 48 49 49 50 50 51 51 52 52 53 53 54 54 55 55 56 56 57 57 58 58 59 59 60 60 61 61 62 62 63 63

しかし,実際の出力はそれぞれ以下の通り.

==================== low ====================
0 0 1 1 2 2 3 3 4 4 5 5 6 6 7 7 16 16 17 17 18 18 19 19 20 20 21 21 22 22 23 23

==================== high ====================
8 8 9 9 10 10 11 11 12 12 13 13 14 14 15 15 24 24 25 25 26 26 27 27 28 28 29 29 30 30 31 31
==================== low ====================
0 0 1 1 2 2 3 3 4 4 5 5 6 6 7 7 16 16 17 17 18 18 19 19 20 20 21 21 22 22 23 23 32 32 33 33 34 34 35 35 36 36 37 37 38 38 39 39 48 48 49 49 50 50 51 51 52 52 53 53 54 54 55 55

==================== high ====================
8 8 9 9 10 10 11 11 12 12 13 13 14 14 15 15 24 24 25 25 26 26 27 27 28 28 29 29 30 30 31 31 40 40 41 41 42 42 43 43 44 44 45 45 46 46 47 47 56 56 57 57 58 58 59 59 60 60 61 61 62 62 63 63

実行結果を観察する限り,AVX,AVX-512のunpackは128bitごとに区切り,SSEのunpack命令を実行しているようになっている.

期待通りの出力を得る,すなわち,レジスタ全体に渡ってインタリーブを行うためには,

  • AVXの場合,unpack命令の後に並び換え
  • AVX-512の場合,8bit毎にレジスタから値を選択

する必要がある.

具体的には以下のようなコードにするとよい. AVXの場合は, _mm256_permute2f128_si256() ,AVX-512の場合は, _mm512_permutex2var_epi8() を用いるとよい. AVXの場合,128bit毎に整列を行う命令が見当たらなかったので,unpackを行わず,直接値を選択する形となる.

// $ g++ -std=gnu++14 -march=native main.cpp -o main
#include <iostream>
#ifdef _MSC_VER
#  include <intrin.h>
#else
#  include <x86intrin.h>
#endif  // _MSC_VER


static inline void
showAVX256(const __m256i& v256) noexcept
{
  alignas(alignof(__m256i)) unsigned char z[sizeof(__m256i)] = {0};
  _mm256_store_si256(reinterpret_cast<__m256i*>(z), v256);
  for (const auto& b : z) {
    std::cout << static_cast<unsigned int>(b) << ' ';
  }
  std::cout << "\n" << std::endl;
}


int
main()
{
  alignas(alignof(__m256i)) unsigned char z[sizeof(__m256i)] = {0};

  for (std::size_t i = 0; i < sizeof(z) / sizeof(z[0]); i++) {
    z[i] = static_cast<unsigned char>(i);
  }

  __m256i v256 = _mm256_load_si256(reinterpret_cast<const __m256i*>(z));
  // unpackの結果を保存
  __m256i vlo = _mm256_unpacklo_epi8(v256, v256);
  __m256i vhi = _mm256_unpackhi_epi8(v256, v256);

  // _mm256_permute2f128_si256() で整列を行う
  std::cout << "==================== permute vlo ====================" << std::endl;
  showAVX256(_mm256_permute2f128_si256(vlo, vhi, 0x20));

  std::cout << "==================== permute vhi ====================" << std::endl;
  showAVX256(_mm256_permute2f128_si256(vlo, vhi, 0x31));

  return 0;
}
// $ g++ -std=gnu++14 -march=native -mavx512vbmi main.cpp -o main
#include <iostream>
#ifdef _MSC_VER
#  include <intrin.h>
#else
#  include <x86intrin.h>
#endif  // _MSC_VER


static inline void
showAVX512(const __m512i& v512) noexcept
{
  alignas(alignof(__m512i)) unsigned char z[sizeof(__m512i)] = {0};
  _mm512_store_si512(reinterpret_cast<__m512i*>(z), v512);
  for (const auto& b : z) {
    std::cout << static_cast<unsigned int>(b) << ' ';
  }
  std::cout << "\n" << std::endl;
}


int
main()
{
  static const __m512i LOIDX = _mm512_setr_epi64(
      0x4303420241014000,
      0x4707460645054404,
      0x4b0b4a0a49094808,
      0x4f0f4e0e4d0d4c0c,
      0x5313521251115010,
      0x5717561655155414,
      0x5b1b5a1a59195818,
      0x5f1f5e1e5d1d5c1c);
  static const __m512i HIIDX = _mm512_setr_epi64(
      0x6323622261216020,
      0x6727662665256424,
      0x6b2b6a2a69296828,
      0x6f2f6e2e6d2d6c2c,
      0x7333723271317030,
      0x7737763675357434,
      0x7b3b7a3a79397838,
      0x7f3f7e3e7d3d7c3c);
  alignas(alignof(__m512i)) unsigned char z[sizeof(__m512i)] = {0};

  for (std::size_t i = 0; i < sizeof(z) / sizeof(z[0]); i++) {
    z[i] = static_cast<unsigned char>(i);
  }

  __m512i v512 = _mm512_load_si512(reinterpret_cast<const __m512i*>(z));

  // _mm512_permutex2var_epi8() で整列を行う
  std::cout << "==================== permute vlo ====================" << std::endl;
  showAVX512(_mm512_permutex2var_epi8(v512, LOIDX, v512));

  std::cout << "==================== permute vhi ====================" << std::endl;
  showAVX512(_mm512_permutex2var_epi8(v512, HIIDX, v512));

  return 0;
}

これで期待通りの出力を得ることが可能である.

_mm256_permute2f128_si256()_mm512_permutex2var_epi8() については,Intelのガイドページより検索するとよい.

をそれぞれ見るとよい. 疑似コードは動作の理解の助けになるはずだ. 特に, _mm256_permute2f128_si256() の第3引数, _mm512_permutex2var_epi8() の第二引数については,疑似コードを見なければ理解しづらい.

_mm256_permute2f128_si256() の第三引数は,下位4byteが出力先の下位4byte,上位4byteが出力先の上位4byteに対応しており,第一引数と第二引数のAVXレジスタを128bitに区切ったセクションのうち,どこから値を取るかを決定する.

なお,AVX-512のインタリーブに関しては,AVXと同様,unpackを行った後,並び換えてもよい. その際,64bit毎に並び替えを行う _mm512_permutex2var_epi64() を用いるとよさそうだ. _mm512_permutex2var_epi8() は8bitごと, _mm512_permutex2var_epi64() は64bitごとに,第一引数と第三引数のAVX-512レジスタのうち,どちらの値を取るかを指定する命令である.

// $ g++ -std=gnu++14 -march=native -mavx512vbmi main.cpp -o main
#include <iostream>
#ifdef _MSC_VER
#  include <intrin.h>
#else
#  include <x86intrin.h>
#endif  // _MSC_VER


static inline void
showAVX512(const __m512i& v512) noexcept
{
  alignas(alignof(__m512i)) unsigned char z[sizeof(__m512i)] = {0};
  _mm512_store_si512(reinterpret_cast<__m512i*>(z), v512);
  for (const auto& b : z) {
    std::cout << static_cast<unsigned int>(b) << ' ';
  }
  std::cout << "\n" << std::endl;
}


int
main()
{
  // 各引数は64bit整数.各引数の2byte目より上位は全て0なので,上位は書いていない.
  static const __m512i LOIDX = _mm512_setr_epi64(0x00, 0x01, 0x08, 0x09, 0x02, 0x03, 0x0a, 0x0b);
  static const __m512i HIIDX = _mm512_setr_epi64(0x04, 0x05, 0x0c, 0x0d, 0x06, 0x07, 0x0e, 0x0f);
  alignas(alignof(__m512i)) unsigned char z[sizeof(__m512i)] = {0};

  for (std::size_t i = 0; i < sizeof(z) / sizeof(z[0]); i++) {
    z[i] = static_cast<unsigned char>(i);
  }

  __m512i v512 = _mm512_load_si512(reinterpret_cast<const __m512i*>(z));
  // unpackの結果を保存
  __m512i vlo = _mm512_unpacklo_epi8(v512, v512);
  __m512i vhi = _mm512_unpackhi_epi8(v512, v512);

  // _mm512_permutex2var_epi64() で整列を行う
  std::cout << "==================== permute vlo ====================" << std::endl;
  showAVX512(_mm512_permutex2var_epi64(vlo, LOIDX, vhi));

  std::cout << "==================== permute vhi ====================" << std::endl;
  showAVX512(_mm512_permutex2var_epi64(vlo, HIIDX, vhi));

  return 0;
}

まとめ

AVX,AVX-512の場合,unpack命令だけではレジスタ全体のインタリーブを行うことができない. unpack後に,並び換えを行う必要がある.

参考文献

SIMDの組み込み関数のことはじめ

C言語 C++

はじめに

現代のCPUではSIMD(Single Instruction Multiple Data)命令を利用することができる. SIMD命令とはその名の通り,ひとつの命令で複数のデータを処理するものである.

Intel系のCPUでは,MMX/SSE/AVX/AVX-512といったSIMD命令が利用可能であり,ARM CPUではNEONというSIMD命令が用意されている. 各SIMDSIMD用のレジスタの対応関係は以下のようになる.

項目 利用可能レジスタ
MMX 64bit のMMレジスタ
SSE 128bit のXMMレジスタ
AVX 256bit のYMMレジシタ
AVX-512 512bit のZMMレジシタ
ARM NEON 64bitのD(Double-Word)レジスタおよび128bitのQ(Quad-Word)レジスタ

これらのレジスタを用いて,例えば4つのint型を一気に処理するといったことを行うのがCPUにおけるSIMDである.

この記事では,このSIMD命令をC/C++から利用することについて記述する.

SIMDをプログラム利用するには

SIMD命令というと小難しそうで,インラインアセンブラを利用しなければならないかというと,そうではない. C/C++から関数の形で利用できるように,各コンパイラで共通のAPIである組み込み関数が提供されている. 組み込み関数とはいえ関数なので,関数呼び出しの形で記述することになるが,実際に関数呼び出しが発生するわけではなく,インライン展開され,対応するアセンブラの命令へとコード生成される.

なお,SIMDレジスタに対して,メモリのロードやストアを行う場合,後述するように利用幅と同じ境界に配置されている位置に対して行う必要がある. 特に,MMX/SSE/AVX/AVX-512の場合,アラインメント条件を満たさなければ,SEGVで落ちる関数がある. 落ちない版の関数もあるが,そういった関数は落ちる関数より動作としては遅い.

ARM NEONは落ちる関数は無いが,アラインメント条件を満たしておいた方が高速に動作すると思われる.

インクルード

何はともあれ,まず組み込み関数が宣言されているヘッダをインクルードしなければ始まらない. 各SIMD命令セットとヘッダの対応関係は以下のようになる.

命令セット ヘッダファイル
MMX <mmintrin.h>
SSE <xmmintrin.h>
SSE2 <emmintrin.h>
SSE3 <pmmintrin.h>
SSSE3 <tmmintrin.h>
SSE4.1 <smmintrin.h>
SSE4.2 <nmmintrin.h>
AES <wmmintrin.h>
AVX, AVX2, FMA <immintrin.h>
AVX-512 <zmmintrin.h>
ARM NEON <arm_neon.h>

MMX/SSE/AVX/AVX-512関連のヘッダは多く,これらをいちいちインクルードするのは面倒である. 現実的にはまとめてインクルードすることが可能なヘッダを利用するのがよい. ただし,MSVCとgcc/clangでヘッダが異なるため,注意しなければならない.

環境 ヘッダファイル
MSVC <intrin.h>
gcc/clang <x86intrin.h>

具体的なインクルード部分のコードを書くと以下のようになる.

#ifdef _MSC_VER
#  include <intrin.h>
#else
#  include <x86intrin.h>
#endif

なお,gcc/clangでも,x64環境ならば <intrin.h> が存在するが,x86環境でも利用可能な方に合わせておく方が何かと都合が良いだろう.

コンパイルオプション

実はヘッダをインクルードするだけではSIMDの組み込み関数は利用できない. 以下のようにコンパイルオプションを指定する必要がある.

命令セット gccのオプション MSVCのオプション 定義されるマクロ
MMX -mmmx /arch:MMX __MMX__
SSE -msse /arch:SSE __SSE__
SSE2 -msse2 /arch:SSE2 __SSE2__
SSE3 -msse3 __SSE3__
SSSE3 -mssse3 __SSSE3__
SSE4.1 -msse4.1 __SSE4_1__
SSE4.2 -msse4.2 __SSE4_2__
AES -maes __AES__
AVX -mavx /arch:AVX __AVX__
AVX2 -mavx2 /arch:AVX2 __AVX2__
FMA -mfma __FMA__
AVX-512 -mavx512* ( *bw, cq, ed など) __AVX512*__
ARM NEON -mfpu=neon など __ARM_NEON または __ARM_NEON__

MMX/SSE/AVX/AVX-512関連のオプションは, -march=native-mtune=native などを指定することで,一括で上記のオプションのうち,利用可能なものを指定できる. ARM CPU環境のgccでは, -march=native-mtune=native と指定することができない場合があり,そのときは利用しているARM CPUに合わせて, -fpu=neon-fp-armv8 などと指定する必要がある(これはRaspberry Pi 3の例).

上記の表では簡略に紹介したが,gccのAVX-512に関するオプションは以下のように多数ある.

  • -mavx512f
  • -mavx512er
  • -mavx512cd
  • -mavx512pf
  • -mavx512dq
  • -mavx512bw
  • -mavx512vl
  • -mavx512ifma
  • -mavx512vbmi

なお,現在のところAVX-512が利用できるCPUは限られている. -march=native を指定したとしても,AVX-512が有効にならない場合の方が多いので,上記のオプションを別途指定すると,コンパイルだけは通るだろう. しかし,非対応のCPUでAVX-512命令を実行したとても,以下のようなエラーメッセージが出力されるだろう. (これはMSYS2でzsh上で実行した結果である)

$ ./main.exe
zsh: illegal hardware instruction  ./main.exe

AVX-512の動作を確認するだけならば,Intel公式のエミュレータを利用するとよい. 予め,AVX-512命令が含まれる実行バイナリを生成し,以下のように実行する.

$ sde -- ./main.exe

変数のアラインメントを指定する

C++11,C11から言語の標準機能として,変数のアラインメントを指定することができるようになったが,それ以前は変数のアラインメントはコンパイラ独自の機能を利用しなければ,指定することができない. 古いコンパイラコンパイルすることを考慮すると,以下のように差を吸収するマクロを定義するとよい.

#include <cstddef>
#include <iostream>

#if defined(__cplusplus) && __cplusplus < 201103L
#  ifdef _MSC_VER
#    define alignas(n)  __declspec(align(n))
#  else
#    define alignas(n)  __attribute__((aligned(n)))
#  endif  // _MSC_VER
#endif  // defined(__cplusplus) && __cplusplus < 201103L


// 以下,利用コード


int
main()
{
  static const int ALIGN = 32;
  alignas(ALIGN) unsigned char array[10] = {0};
  if ((reinterpret_cast<std::ptrdiff_t>(array)) % ALIGN == 0) {
    std::cout << "Static array is " << ALIGN << " byte aligned.\n";
  } else {
    std::cout << "Static array is not " << ALIGN << " byte aligned.\n";
  }

  return 0;
}

上記はC++用だが,C言語なら以下のように定義するとよい.

#include <stddef.h>
#include <stdio.h>

#if defined(__STDC_VERSION__) && __STDC_VERSION__ < 201102L
#  ifdef _MSC_VER
#    define _Alignas(n)  __declspec(align(n))
#  else
#    define _Alignas(n)  __attribute__((aligned(n)))
#  endif  // _MSC_VER
#endif  // defined(__cplusplus) && __cplusplus < 201103L


/* 以下,利用コード */


#define ALIGN  32

int
main(void)
{
  _Alignas(ALIGN) unsigned char array[10] = {0};
  if ((ptrdiff_t) array % ALIGN == 0) {
    printf("Static array is %d byte aligned.\n", ALIGN);
  } else {
    printf("Static array is not %d byte aligned.\n", ALIGN);
  }

  return 0;
}

アラインされたメモリを動的確保する

通常のC/C++における std::malloc()std::calloc()new 等では16byteや32byte境界にアラインメントされたメモリを動的確保することはできない. 以下に示す専用のメモリ確保関数が必要となる. (C言語の場合, <cstdlib><stdlib.h> に読み換えること)

メモリ確保関数 メモリ解放関数 ヘッダ 特徴
_aligned_malloc() _aligned_free() <malloc.h> MSVCのみ.
posix_memalign() std::free() <cstdlib> gcc/clangのみ.
aligned_alloc() std::free() <cstdlib> gcc/clangのみ.確保サイズはアラインメントの整数倍に限る
memalign() std::free() <malloc.h> gcc/clangのみ.廃止されているとのこと.
_mm_malloc() _mm_free() <malloc.h> Intel CPUのみ.

種々のアラインされたメモリ確保関数があり,どれを利用すればいいか判断に困るかもしれない. しかし,おおまかには,以下のように利用する関数を判断すればよい.

  • MSVCなら _aligned_malloc()_aligned_free()
  • gcc/clangなら posix_memalign()std::free()

これを考慮し,条件コンパイルで利用する関数を分岐するラッパー関数を作るとよい. 簡単なコードは以下のようになる.

なお,C++11以降, std::align()std::aligned_storage() といった関数が利用できるが, std::align() は既に確保されたバッファの指定されたアドレスからポインタを進め,アラインメント条件を満たす位置のアドレスを返却するだけの関数であり, std::aligned_storage() はアラインされた静的配列を作成するための関数なので,やや使い勝手が悪いといえる.

// <type_traits> はC++11以降のものなので,それ以前でコンパイルしたい場合は関連部分を削除すること
#include <cstddef>
#include <iostream>
#include <memory>
#include <type_traits>
#if defined(_MSC_VER) || defined(__MINGW32__)
#  include <malloc.h>
#else
#  include <cstdlib>
#endif  // defined(_MSC_VER) || defined(__MINGW32__)


/*!
 * @brief アラインメントされたメモリを動的確保する関数
 * @param [in] size       確保するメモリサイズ (単位はbyte)
 * @param [in] alignment  アラインメント (2のべき乗を指定すること)
 * @return  アラインメントし,動的確保されたメモリ領域へのポインタ
 */
template<typename T = void*, typename std::enable_if<std::is_pointer<T>::value, std::nullptr_t>::type = nullptr>
static inline T
alignedMalloc(std::size_t size, std::size_t alignment) noexcept
{
#if defined(_MSC_VER) || defined(__MINGW32__)
  return reinterpret_cast<T>(_aligned_malloc(size, alignment));
#else
  void* p;
  return reinterpret_cast<T>(posix_memalign(&p, alignment, size) == 0 ? p : nullptr);
#endif  // defined(_MSC_VER) || defined(__MINGW32__)
}


/*!
 * @brief アラインメントされたメモリを解放する関数
 * @param [in] ptr  解放対象のメモリの先頭番地を指すポインタ
 */
static inline void
alignedFree(void* ptr) noexcept
{
#if defined(_MSC_VER) || defined(__MINGW32__)
  _aligned_free(ptr);
#else
  std::free(ptr);
#endif  // defined(_MSC_VER) || defined(__MINGW32__)
}


// 以下,利用コード


/*!
 * @brief std::unique_ptr で利用するアラインされたメモリ用のカスタムデリータ
 */
struct AlignedDeleter
{
  void
  operator()(void* p) const noexcept
  {
    alignedFree(p);
  }
};


int
main()
{
  static constexpr int ALIGN = 32;
  std::unique_ptr<unsigned char[], AlignedDeleter> array(alignedMalloc<unsigned char*>(10 * sizeof(unsigned char), ALIGN));
  if (array.get() == nullptr) {
    std::cerr << "Failed to allocate memory" << std::endl;
    return 1;
  }

  if ((reinterpret_cast<std::ptrdiff_t>(array.get())) % ALIGN == 0) {
    std::cout << "Dynamic allocated memory is " << ALIGN << " byte aligned.\n";
  } else {
    std::cout << "Dynamic allocated memory is not " << ALIGN << " byte aligned.\n";
  }

  return 0;
}

このコードはC++11の範疇のものであるが,C言語の範囲で書き直すと以下のようになる. C99以降は inline が利用可能であるが,古いコンパイラを使用することを考慮し,置き換えるマクロを記述する.

#include <stdio.h>
#include <stddef.h>
#if defined(_MSC_VER) || defined(__MINGW32__)
#  include <malloc.h>
#else
#  include <stdlib.h>
#endif  /* defined(_MSC_VER) || defined(__MINGW32__) */

#ifndef __cplusplus
#  if defined(_MSC_VER)
#    define inline      __inline
#    define __inline__  __inline
#  elif !defined(__GNUC__) && !defined(__STDC_VERSION__) || __STDC_VERSION__ < 199901L
#    define inline
#    define __inline
#  endif
#endif


/*!
 * @brief アラインメントされたメモリを動的確保する関数
 * @param [in] size       確保するメモリサイズ (単位はbyte)
 * @param [in] alignment  アラインメント (2のべき乗を指定すること)
 * @return  アラインメントし,動的確保されたメモリ領域へのポインタ
 */
static inline void*
alignedMalloc(size_t size, size_t alignment)
{
#if defined(_MSC_VER) || defined(__MINGW32__)
  return _aligned_malloc(size, alignment);
#else
  void* p;
  return posix_memalign((void**) &p, alignment, size) == 0 ? p : NULL;
#endif  /* _MSC_VER */
}


/*!
 * @brief アラインメントされたメモリを解放する関数
 * @param [in] ptr  解放対象のメモリの先頭番地を指すポインタ
 */
static inline void
alignedFree(void* ptr)
{
#if defined(_MSC_VER) || defined(__MINGW32__)
  _aligned_free(ptr);
#else
  free(ptr);
#endif  /* _MSC_VER */
}


/* 以下,利用コード */


int
main(void)
{
  static const int ALIGN = 32;
  unsigned char* array = (unsigned char*) alignedMalloc(10 * sizeof(unsigned char), ALIGN);
  if (array == NULL) {
    fprintf(stderr, "Failed to allocate memory\n");
    return 1;
  }

  if (((ptrdiff_t) array) % ALIGN == 0) {
    printf("Dynamic allocated memory is %d byte aligned.\n", ALIGN);
  } else {
    printf("Dynamic allocated memory is not %d byte aligned.\n", ALIGN);
  }
  alignedFree(array);

  return 0;
}

なお,今回は適当なアラインメントを指定したが,実際にSSE/AVX/AVX-512/NEONを用いるときは,SSE/AVX/AVX-512/NEONの変数型からアラインメントを取得するとよい. 型や変数からアラインメントを得る機能はC++11およびC11以降であれば, alignof 演算子で取得でき,それ以前の環境であれば,コンパイラ拡張機能を用いることで取得できる. この差を吸収するなら,以下のようなマクロを定義するとよい.

#if defined(__cplusplus) && __cplusplus < 201103L
#  ifdef _MSC_VER
#    define alignof(n)  __alignof(n)
#  else
#    define alignof(n)  __alignof__(n)
#  endif  // _MSC_VER
#endif  // defined(__cplusplus) || cplusplus < 201103L

// alingas(alignof(__m256i)) のような形で使用

SSE/AVX/NEON のサンプルコード

簡単なサンプルコードをSSE/AVX/NEONの例として提示する. このコードはMSVC/gc/clangのいずれのコンパイラでもコンパイルすることができるようにしている.

AVX-512については,利用可能なCPUを搭載したマシンが手元に無いため割愛するが,AVXと同様のコードで記述できると思う. コンパイル時に以下のマクロを定義すると,対応した命令を用いたコードが有効化される. 有効化しようとしても,コンパイラが対応していない場合は,冒頭の部分でエラーが発生するはずだ. また,以下のいずれのマクロも定義しなかった場合,SIMDを用いないコードとなる.

マクロ 有効化されるSIMD
ENABLE_AVX AVX
ENABLE_SSE SSE
ENABLE_NEON ARM NEON

直接的に __AVX____SSE2__ 等のマクロが定義されているかどうかで判断しないのは,コンパイラAPIとして提供していたとしても,CPUが対応しておらず,SIMD命令を利用できない場合もあるからだ. また,AVXやSSEの切り替えが容易になり,ベンチマークテストがしやすいという利点もあるだろう.

さて,具体的には以下のようにオプションを指定してコンパイルするとよい. gccの場合は,

有効化する機能 コマンド
AVX-512 $ g++ -std=gnu++11 -march=native -mavx512f -DENABLE_AVX main.cpp -o main.o
AVX $ g++ -std=gnu++11 -march=native -DENABLE_AVX main.cpp -o main.o
SSE $ g++ -std=gnu++11 -march=native -DENABLE_SSE main.cpp -o main.o
ARM NEON $ g++ -std=gnu++11 -mfpu=neon-fp-armv8 -DENABLE_NEON main.cpp -o main.o
SIMDを利用しない $ g++ -std=gnu++11 main.cpp -o main.o

であり,MSVCの場合は,

有効化する機能 コマンド
AVX > cl.exe /arch:AVX /DENABLE_AVX main.cpp
SSE > cl.exe /arch:SSE2 /DENABLE_SSE main.cpp
SIMDを利用しない > cl.exe main.cpp

といった具合である.

ベクトルの内積計算

定番のベクトルの内積を計算するコードを示す. FMA(積和演算)が利用可能な場合は,そちらを用いて,高速に処理できるようにしてある. また,SIMDを用いない場合であっても,C++11/C11以降で <cmath> から提供されている std::fma() を用いることで,内積計算の高速化が期待できるようにする.

#if defined(ENABLE_AVX512) && !defined(__AVX512F__)
#  error Macro: ENABLE_AVX512 is defined, but unable to use AVX512F intrinsic functions
#elif defined(ENABLE_AVX) && !defined(__AVX__)
#  error Macro: ENABLE_AVX is defined, but unable to use AVX intrinsic functions
#elif defined(ENABLE_SSE) && !defined(__SSE2__)
#  error Macro: ENABLE_SSE is defined, but unable to use SSE intrinsic functions
#elif defined(ENABLE_NEON) && !defined(__ARM_NEON) && !defined(__ARM_NEON__)
#  error Macro: ENABLE_NEON is defined, but unable to use NEON intrinsic functions
#else


#include <cmath>
#include <cstddef>
#include <iostream>
#include <memory>
#include <type_traits>
#if defined(_MSC_VER) || defined(__MINGW32__)
#  include <malloc.h>
#else
#  include <cstdlib>
#endif
#if defined(ENABLE_AVX512) || defined(ENABLE_AVX) || defined(ENABLE_SSE)
#  ifdef _MSC_VER
#    include <intrin.h>
#  else
#    include <x86intrin.h>
#  endif  // _MSC_VER
#elif defined(ENABLE_NEON)
#  include <arm_neon.h>
#endif  // defined(ENABLE_AVX512) || defined(ENABLE_AVX) || defined(ENABLE_SSE)

#if defined(__cplusplus) && __cplusplus < 201103L
#  ifdef _MSC_VER
#    define alignof(n)  __alignof(n)
#  else
#    define alignof(n)  __alignof(n)__
#  endif  // _MSC_VER
#endif  // defined(__cplusplus) || cplusplus < 201103L


/*!
 * @brief アラインメントされたメモリを動的確保する関数
 * @param [in] size       確保するメモリサイズ (単位はbyte)
 * @param [in] alignment  アラインメント (2のべき乗を指定すること)
 * @return  アラインメントし,動的確保されたメモリ領域へのポインタ
 */
template<typename T = void*, typename std::enable_if<std::is_pointer<T>::value, std::nullptr_t>::type = nullptr>
static inline T
alignedMalloc(std::size_t size, std::size_t alignment) noexcept
{
#if defined(_MSC_VER) || defined(__MINGW32__)
  return reinterpret_cast<T>(_aligned_malloc(size, alignment));
#else
  void* p;
  return reinterpret_cast<T>(posix_memalign(&p, alignment, size) == 0 ? p : nullptr);
#endif  // defined(_MSC_VER) || defined(__MINGW32__)
}


/*!
 * @brief アラインメントされたメモリを解放する関数
 * @param [in] ptr  解放対象のメモリの先頭番地を指すポインタ
 */
static inline void
alignedFree(void* ptr) noexcept
{
#if defined(_MSC_VER) || defined(__MINGW32__)
  _aligned_free(ptr);
#else
  std::free(ptr);
#endif  // defined(_MSC_VER) || defined(__MINGW32__)
}


/*!
 * @brief std::unique_ptr で利用するアラインされたメモリ用のカスタムデリータ
 */
struct AlignedDeleter
{
  void
  operator()(void* p) const noexcept
  {
    alignedFree(p);
  }
};


#if defined(ENABLE_AVX512)
static constexpr int ALIGN = alignof(__m512);
#elif defined(ENABLE_AVX)
static constexpr int ALIGN = alignof(__m256);
#elif defined(ENABLE_SSE)
static constexpr int ALIGN = alignof(__m128);
#elif defined(ENABLE_NEON)
static constexpr int ALIGN = alignof(float32x4_t);
#else
static constexpr int ALIGN = 8;
#endif  // defined(ENABLE_AVX512)


/*!
 * @brief 内積計算を行う関数
 * @param [in] a  ベクトルその1
 * @param [in] b  ベクトルその2
 * @param [in] n  ベクトルのサイズ
 * @return  内積
 */
static inline float
innerProduct(const float* a, const float* b, std::size_t n)
{
#if defined(ENABLE_AVX512)
  static constexpr std::size_t INTERVAL = sizeof(__m512) / sizeof(float);
  __m512 sumx16 = {0};
  for (std::size_t i = 0; i < n; i += INTERVAL) {
    __m512 ax16 = _mm512_load_ps(&a[i]);
    __m512 bx16 = _mm512_load_ps(&b[i]);
#  ifdef __FMA__
    sumx16 = _mm512_fmadd_ps(ax16, bx16, sumx16);
#  else
    ax16 = _mm512_mul_ps(ax16, bx16);
    sumx16 = _mm512_add_ps(sumx16, ax16);
#  endif  // __FMA__
  }

  alignas(ALIGN) float s[INTERVAL] = {0};
  _mm512_store_ps(s, sumx16);
  float sum = s[0] + s[1] + s[2] + s[3] + s[4] + s[5] + s[6] + s[7]
    + s[8] + s[9] + s[10] + s[11] + s[12] + s[13] + s[14] + s[15];

  for (std::size_t i = n - n % INTERVAL; i < n; i++) {
    sum += a[i] * b[i];
  }
  return sum;
#elif defined(ENABLE_AVX)
  static constexpr std::size_t INTERVAL = sizeof(__m256) / sizeof(float);
  __m256 sumx8 = {0};
  for (std::size_t i = 0; i < n; i += INTERVAL) {
    __m256 ax8 = _mm256_load_ps(&a[i]);
    __m256 bx8 = _mm256_load_ps(&b[i]);
#  ifdef __FMA__
    sumx8 = _mm256_fmadd_ps(ax8, bx8, sumx8);
#  else
    ax8 = _mm256_mul_ps(ax8, bx8);
    sumx8 = _mm256_add_ps(sumx8, ax8);
#  endif  // __FMA__
  }

  alignas(ALIGN) float s[INTERVAL] = {0};
  _mm256_store_ps(s, sumx8);
  float sum = s[0] + s[1] + s[2] + s[3] + s[4] + s[5] + s[6] + s[7];

  for (std::size_t i = n - n % INTERVAL; i < n; i++) {
    sum += a[i] * b[i];
  }
  return sum;
#elif defined(ENABLE_SSE)
  static constexpr std::size_t INTERVAL = sizeof(__m128) / sizeof(float);
  __m128 sumx4 = {0};
  for (std::size_t i = 0; i < n; i += INTERVAL) {
    __m128 ax4 = _mm_load_ps(&a[i]);
    __m128 bx4 = _mm_load_ps(&b[i]);
#  ifdef __FMA__
    sumx4 = _mm_fmadd_ps(ax4, bx4, sumx4);
#  else
    ax4 = _mm_mul_ps(ax4, bx4);
    sumx4 = _mm_add_ps(sumx4, ax4);
#  endif  // __FMA__
  }

  alignas(ALIGN) float s[INTERVAL] = {0};
  _mm_store_ps(s, sumx4);
  float sum = s[0] + s[1] + s[2] + s[3];

  for (std::size_t i = n - n % INTERVAL; i < n; i++) {
    sum += a[i] * b[i];
  }
  return sum;
#elif defined(ENABLE_NEON)
  static constexpr std::size_t INTERVAL = sizeof(float32x4_t) / sizeof(float);
  float32x4_t sumx4 = {0};
  for (std::size_t i = 0; i < n; i += INTERVAL) {
    float32x4_t ax4 = vld1q_f32(&a[i]);
    float32x4_t bx4 = vld1q_f32(&b[i]);
    sumx4 = vmlaq_f32(sumx4, ax4, bx4);
  }

  float sum = sumx.val[0] + sumx.val[1] + sumx.val[2] + sumx.val[3];

  for (std::size_t i = n - n % INTERVAL; i < n; i++) {
    sum += a[i] * b[i];
  }
  return sum;
#else
  float sum = 0.0;
  for (std::size_t i = 0; i < n; i++) {
    // <cmath>のstd::fma関数を用いると,積和演算がハードウェアのサポートを受けることを期待できる
    // 処理としては, sum += a[i] * b[i]; と同じ
    sum = std::fma(a[i], b[i], sum);
  }
  return sum;
#endif  // defined(ENABLE_AVX512)
}


int
main()
{
  static constexpr int N_ELEMENT = 256;

  std::unique_ptr<float[], AlignedDeleter> a(alignedMalloc<float*>(N_ELEMENT * sizeof(float), ALIGN));
  std::unique_ptr<float[], AlignedDeleter> b(alignedMalloc<float*>(N_ELEMENT * sizeof(float), ALIGN));
  for (int i = 0; i < N_ELEMENT; i++) {
    a[i] = static_cast<float>(i);
    b[i] = static_cast<float>(i);
  }
  std::cout << innerProduct(a.get(), b.get(), N_ELEMENT) << std::endl;

  return 0;
}


#endif  // defined(ENABLE_AVX512) && !defined(__AVX512F__)

最近傍法による画像の2倍拡大

最近傍法,すなわち単純なピクセルコピーのみを行って,8bitグレースケール画像を2倍に拡大するコードを記述する. 2倍拡大という条件に限定すれば,出力先画像のインデックス値のとる値が単純になるので,SIMDで簡単に処理を記述できる.

読み込む画像ファイル名は test.jpg とし,読み込みにOpenCVを用いる. 画像ファイルの横幅は,16または32の倍数でなければならない.

コンパイルは以下のようにするとよい.

$ g++ -std=gnu++11 main.cpp -march=native -DENABLE_AVX -I/usr/include/opencv -I/usr/include/opencv2 -lopencv_core -lopencv_highgui -lopencv_imgcodecs -o main.o

AVX-512を利用する場合は, -mavx512vbmi -DENABLE_AVX512 を付加するとよい.

なお,OpenCVcv::Matカスタムアロケータを適用することができるらしいが,コードが煩雑になりそうなので,SSE/AVXにおいてはアラインメント条件を満たさなくてもよい関数を用いている.

#if defined(ENABLE_AVX512) && !defined(__AVX512F__)
#  error Macro: ENABLE_AVX512 is defined, but unable to use AVX512F intrinsic functions
#elif defined(ENABLE_AVX) && !defined(__AVX__)
#  error Macro: ENABLE_AVX is defined, but unable to use AVX intrinsic functions
#elif defined(ENABLE_SSE) && !defined(__SSE2__)
#  error Macro: ENABLE_SSE is defined, but unable to use SSE intrinsic functions
#elif defined(ENABLE_NEON) && !defined(__ARM_NEON) && !defined(__ARM_NEON__)
#  error Macro: ENABLE_NEON is defined, but unable to use NEON intrinsic functions
#else  // defined(ENABLE_AVX512) && !defined(__AVX512F__)


#include <cmath>
#include <cstddef>
#include <iostream>
#include <memory>
#include <type_traits>
#if defined(_MSC_VER) || defined(__MINGW32__)
#  include <malloc.h>
#else
#  include <cstdlib>
#endif
#if defined(ENABLE_AVX512) || defined(ENABLE_AVX) || defined(ENABLE_SSE)
#  ifdef _MSC_VER
#    include <intrin.h>
#  else
#    include <x86intrin.h>
#  endif  // _MSC_VER
#elif defined(ENABLE_NEON)
#  include <arm_neon.h>
#endif  // defined(ENABLE_AVX512) || defined(ENABLE_AVX) || defined(ENABLE_SSE)

#include <opencv2/opencv.hpp>


#if defined(_MSC_VER) && _MSC_VER >= 1400 || \
  defined(__GNUC__) && defined(__GNUC_MINOR__) && (__GNUC__ > 2 || __GNUC__ == 2 && __GNUC_MINOR__ >= 92)
#  define restrict  __restrict
#else
#  define restrict
#endif


/*!
 * @brief アラインメントされたメモリを動的確保する関数
 * @param [in] size       確保するメモリサイズ (単位はbyte)
 * @param [in] alignment  アラインメント (2のべき乗を指定すること)
 * @return  アラインメントし,動的確保されたメモリ領域へのポインタ
 */
template<typename T = void*, typename std::enable_if<std::is_pointer<T>::value, std::nullptr_t>::type = nullptr>
static inline T
alignedMalloc(std::size_t size, std::size_t alignment) noexcept
{
#if defined(_MSC_VER) || defined(__MINGW32__)
  return reinterpret_cast<T>(_aligned_malloc(size, alignment));
#else
  void* p;
  return reinterpret_cast<T>(posix_memalign(&p, alignment, size) == 0 ? p : nullptr);
#endif  // defined(_MSC_VER) || defined(__MINGW32__)
}


/*!
 * @brief アラインメントされたメモリを解放する関数
 * @param [in] ptr  解放対象のメモリの先頭番地を指すポインタ
 */
static inline void
alignedFree(void* ptr) noexcept
{
#if defined(_MSC_VER) || defined(__MINGW32__)
  _aligned_free(ptr);
#else
  std::free(ptr);
#endif  // defined(_MSC_VER) || defined(__MINGW32__)
}


#if defined(ENABLE_AVX512)
static constexpr int ALIGN = alignof(__m512i);
#elif defined(ENABLE_AVX)
static constexpr int ALIGN = alignof(__m256i);
#elif defined(ENABLE_SSE)
static constexpr int ALIGN = alignof(__m128i);
#elif defined(ENABLE_NEON)
static constexpr int ALIGN = alignof(uint8x16_t);
#else
static constexpr int ALIGN = 8;
#endif  // defined(ENABLE_AVX512)


/*!
 * @brief 入力画像データを最近傍法により,2倍のサイズに拡大する
 * @param [out] dstImageData  出力画像データ領域の先頭へのポインタ
 * @param [in]  dstWidth      出力画像データの横幅
 * @param [in]  dstHeight     出力画像データの縦幅
 * @param [in]  srcImageData  入力画像データ領域の先頭へのポインタ
 * @param [in]  srcWidth      入力画像データの横幅
 * @param [in]  srcHeight     入力画像データの縦幅
 * @return  アラインメントし,動的確保されたメモリ領域へのポインタ
 */
static inline void
scale2x(
    unsigned char* restrict dstImageData,
    int dstWidth,
    int dstHeight,
    const unsigned char* restrict srcImageData,
    int srcWidth,
    int srcHeight) noexcept
{
  static constexpr int X_RATIO = 2;
  static constexpr int Y_RATIO = 2;
#if defined(ENABLE_AVX512)
  static constexpr int INTERVAL = sizeof(__m512i) / sizeof(unsigned char);
  static const __m512i LOWIDX = _mm512_setr_epi64(
      0x4303420241014000,
      0x4707460645054404,
      0x4b0b4a0a49094808,
      0x4f0f4e0e4d0d4c0c,
      0x5313521251115010,
      0x5717561655155414,
      0x5b1b5a1a59195818,
      0x5f1f5e1e5d1d5c1c);
  static const __m512i HIGHIDX = _mm512_setr_epi64(
      0x6323622261216020,
      0x6727662665256424,
      0x6b2b6a2a69296828,
      0x6f2f6e2e6d2d6c2c,
      0x7333723271317030,
      0x7737763675357434,
      0x7b3b7a3a79397838,
      0x7f3f7e3e7d3d7c3c);
#elif defined(ENABLE_AVX)
  static constexpr int INTERVAL = sizeof(__m256i) / sizeof(unsigned char);
#elif defined(ENABLE_SSE)
  static constexpr int INTERVAL = sizeof(__m128i) / sizeof(unsigned char);
#elif defined(ENABLE_NEON)
  static constexpr int INTERVAL = sizeof(uint8x16_t) / sizeof(unsigned char);
#else
  static constexpr int INTERVAL = sizeof(unsigned char);
#endif  // defined(ENABLE_AVX512)

  for (int i = 0; i < dstHeight; i++) {
    for (int j = 0; j < dstWidth; j += INTERVAL * X_RATIO) {
#if defined(ENABLE_AVX512)
      // 64pixel分の画素データをロード
      __m512i v512 = _mm512_loadu_si512(reinterpret_cast<const __m512i*>(&srcImageData[i / Y_RATIO * srcWidth + j / X_RATIO]));
      // インタリーブ
      __m512i v512l = _mm512_permutex2var_epi8(v512, LOWIDX, v512);
      __m512i v512u = _mm512_permutex2var_epi8(v512, HIGHIDX, v512);
      // 64pixel x 2のデータを書き込み
      _mm512_storeu_si512(reinterpret_cast<__m512i*>(&dstImageData[i * dstWidth + j + sizeof(__m512i) * 0]), v512l);
      _mm512_storeu_si512(reinterpret_cast<__m512i*>(&dstImageData[i * dstWidth + j + sizeof(__m512i) * 1]), v512u);
#elif defined(ENABLE_AVX)
      // 32pixel分の画素データをロード
      __m256i v256 = _mm256_loadu_si256(reinterpret_cast<const __m256i*>(&srcImageData[i / Y_RATIO * srcWidth + j / X_RATIO]));
      // インタリーブ
      __m256i v256l_ = _mm256_unpacklo_epi8(v256, v256);
      __m256i v256u_ = _mm256_unpackhi_epi8(v256, v256);
      // 上下128bit交換
      __m256i v256l = _mm256_permute2f128_si256(v256l_, v256u_, 0x20);
      __m256i v256u = _mm256_permute2f128_si256(v256l_, v256u_, 0x31);
      // 32pixel x 2のデータを書き込み
      _mm256_storeu_si256(reinterpret_cast<__m256i*>(&dstImageData[i * dstWidth + j + sizeof(__m256i) * 0]), v256l);
      _mm256_storeu_si256(reinterpret_cast<__m256i*>(&dstImageData[i * dstWidth + j + sizeof(__m256i) * 1]), v256u);
#elif defined(ENABLE_SSE)
      // 16pixel分の画素データをロード
      __m128i v128 = _mm_loadu_si128(reinterpret_cast<const __m128i*>(&srcImageData[i / Y_RATIO * srcWidth + j / X_RATIO]));
      // インタリーブ
      __m128i v128l = _mm_unpacklo_epi8(v128, v128);
      __m128i v128u = _mm_unpackhi_epi8(v128, v128);
      // 16pixel x 2のデータを書き込み
      _mm_storeu_si128(reinterpret_cast<__m128i*>(&dstImageData[i * dstWidth + j + sizeof(__m128i) * 0]), v128l);
      _mm_storeu_si128(reinterpret_cast<__m128i*>(&dstImageData[i * dstWidth + j + sizeof(__m128i) * 1]), v128u);
#elif defined(ENABLE_NEON)
      // 16pixel分の画素データをロード
      uint8x16_t v128 = vld1q_u8(&srcImageData[i / Y_RATIO * srcWidth + j]);
      // インタリーブ
      uint8x16x2_t v128x2 = vzipq_u8(v128, v128);
      // 16pixel x 2のデータを書き込み
      vst1q_u8(dstImageData[i * dstWidth + j + sizeof(uint8x16_t) * 0], v128x2.val[0]);
      vst1q_u8(dstImageData[i * dstWidth + j + sizeof(uint8x16_t) * 1], v128x2.val[1]);
#else
      dstImageData[i * dstWidth + j] = srcImageData[i / Y_RATIO * srcWidth + j];
#endif  // defined(ENABLE_AVX512)
    }
  }
}


int
main()
{
  cv::Mat img = cv::imread("test.jpg", 0);
  if (img.data == nullptr) {
    std::cerr << "Cannot open image file: test.jpg" << std::endl;
    return 1;
  }
  cv::Mat scaledImg(cv::Size(img.cols * 2, img.rows * 2), CV_8UC1);
  scale2x(scaledImg.data, scaledImg.cols, scaledImg.rows, img.data, img.cols, img.rows);

  cv::namedWindow("src", CV_WINDOW_AUTOSIZE);
  cv::namedWindow("scaled", CV_WINDOW_AUTOSIZE);
  cv::imshow("src", img);
  cv::imshow("scaled", scaledImg);
  std::cout << "Please hit any key on the window to exit this program" << std::endl;
  cv::waitKey(0);

  return 0;
}


#endif  // defined(ENABLE_AVX512) && !defined(__AVX512F__)

何の脈略も無しに,SSE/AVXやARM NEONの組み込み関数や型を利用したが,SSE/AVXに関してはIntelIntrinsics Guideを,ARM NEONに関してはARM NEON Intrinsicsを参照するとよい.

SSE/AVXの変数型は以下の通り.

内容
__m128 float 型4個分
__m128d double 型2個分
__m128i 整数型 (intunsigned char などを格納できる)
__m256 float 型8個分
__m256d double 型4個分
__m256i 整数型 (intunsigned char などを格納できる)
__m512 float 型16個分
__m512d double 型8個分
__m512i 整数型 (intunsigned char などを格納できる)

SSE/AVXの組み込み関数は基本的に

  • SSEの場合, _mm_[xxx]{[u]}_[yyy]
  • AVXの場合, _mm256_[xxx]{[u]}_[yyy]
  • AVX-512の場合, _mm512_[xxx]{[u]}_[yyy]

の形式で命名されている. [xxx][{u}][yyy] の部分については以下の通り.

該当部分 内容
[xxx] loadstore など,行いたい命令がここにくる
[u] u が付いている関数はアラインメント条件を満たしていなくても,SEGVで落ちない
[yyy] 引数の型によって変化する. ps なら __m128pd なら __m128dsi128 なら __m128i

pspd はそれぞれ Precision Single, Precision Double の略であるそうだ. (si は調べていない)

ARM NEONの変数型は見た目通り, [xxx][size]x[NNN]{x[MMM]} の形式となっている.

該当部分 内容
[xxx] uintintfloat などのベクタの1要素の型がここにくる
[size] ベクタの要素型1つのサイズ (単位はbit)
[NNN] ベクタ要素の個数
[MMM] インタリーブ用にくっつけたNEONレジスタの個数.2から4までの値ろ取り,1つの場合は省略される

ARM NEONの組み込み関数も直感的に利用できる命名で, v[xxx]{[q]}_{yyy} となっている.

該当部分 内容
[xxx] addld など,行いたい命令がここにくる
[q] qが付いていればQレジスタ(128bit)を用いる命令,付いていないならばDレジスタ(64bit)を用いる命令
[yyy] 引数の型によって変化する. u8s16f32 など

まとめ

この記事では以下のことを紹介した.

  • SIMDの概要
  • SIMDの組み込み関数の利用方法
  • コンパイラの差を吸収するアラインメントの指定方法
  • ベクトルの内積を計算するサンプルコード

特に,SIMDの組み込み関数の利用方法を簡単にまとめると以下のようになる.

  • alignas(alignof(__m256i)) ... の形で,変数のアラインメント指定
    • 古いMSVCなら __declspec(align(32))
    • 古いgccなら __attribute__((aligned(32)))
  • gcc
    • #include <x86intrin>
    • $ g++ -march=native ...
    • pisix_memalign() でアラインされた動的メモリ確保, std::free() で解放
  • MSVC
    • #include <intrin>
    • > cl.exe /arch:AVX2 ...
    • _aligned_malloc() でアラインされた動的メモリ確保, _aligned_free() で解放
  • AVX-512非対応のCPUでAVX-512をテストする場合は,Intelのエミュレータを利用

この記事はあくまでSIMDの基礎に過ぎないが,あとは組み込み関数を調べ,うまく組み合わせることで,SIMDをプログラムに組み込めるようになるかもしれない.

参考文献