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

koturnの日記

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

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をプログラムに組み込めるようになるかもしれない.

参考文献

GNU MakeでC/C++の依存関係を自動定義する

背景

C/C++用のMakefileを書いていると,依存関係を記述するのがなかなか厄介に感じられる. 何とか自動的に依存関係を定義したいと考えた.

gcc/g++による依存関係の抽出

gcc/g++は以下のようにして依存関係を抽出することが可能である.

$ g++ -MM foo.cpp
foo.o: foo.cpp aaa.h bbb.h

この機能を利用して,Makefileで依存関係の自動定義を行う.

外部ファイルに依存関係を書き出す

makeには他のMakefileをインクルードすることが可能である. これを利用し,依存関係を記述したファイルを掃き出し,Makefileからインクルードするとよさそうだ.

CXX     := g++
SRCS    := foo.cpp bar.cpp
DEPENDS := depends.mk

.PHONY: depends

depends:
    $(CXX) -MM $(SRCS) > $(DEPENDS)

-include $(DEPENDS)

ただし,この方法だと,make depends を実行しなければ,依存関係を更新されない. また,初回はdepends.mkが存在せず,依存関係が定義されない. これを,自動的に更新したい.

外部ファイルを用いず,依存関係を自動定義する

GNU Makeには,引数をシェルコマンドとして実行し,出力結果として返却する shell 関数,引数をMakefileに記述されたものとして扱う eval 関数が存在する. これらを利用することで,Makefileから依存関係を定義することができそうだ.

CXX  := g++
SRCS := foo.cpp bar.cpp

$(eval $(shell $(CXX) -MM $(SRCS)))

しかし,この方法では改行が無視されるため,うまくいかない. 具体的には,複数ファイルを一度に g++ -MM の引数として渡せない,行継続の \ によってエラーとなる. したがって,ソースファイルを1つずつに対し, g++ -MM を実行する必要し, \ を消去する必要がある.

GNU Makeには,文字列の置換を行う subst 関数, 与えられたリストの要素それぞれに,指定した処理を実行する foreach 関数がある. これらを利用し,前述の問題を解決する.

CXX  := g++
SRCS := foo.cpp bar.cpp

$(foreach SRC,$(SRCS),$(eval $(subst \,,$(shell $(CXX) -MM $(SRC)))))

具体的なMakefileの全貌は以下のようになるだろう. 圧倒的に依存関係の記述が楽になり,ヘッダファイル名は全く記述する必要が無くなっていることがわかる.

CXX      := g++ -std=gnu++14
CXXFLAGS := -Wall -Wextra -O3 -march=native -DNDEBUG
LDFLAGS  := -s
TARGET   := main
OBJS     := $(addsuffix .o, $(basename $(TARGET)) foo bar)
SRCS     := $(OBJS:.o=.cpp)

ifeq ($(OS),Windows_NT)
    TARGET := $(addsuffix .exe, $(TARGET))
else
    TARGET := $(addsuffix .out, $(TARGET))
endif

%.exe:
  $(CXX) $(LDFLAGS) $(filter %.c %.cpp %.cxx %.cc %.o, $^) $(LDLIBS) -o $@
%.out:
  $(CXX) $(LDFLAGS) $(filter %.c %.cpp %.cxx %.cc %.o, $^) $(LDLIBS) -o $@

.PHONY: all clean

all: $(TARGET)
$(TARGET): $(OBJS)

$(foreach SRC,$(SRCS),$(eval $(subst \,,$(shell $(CXX) -MM $(SRC)))))

clean:
    $(RM) $(TARGET) $(OBJS)

まとめ

GNU Make用のMakefileに以下のような記述を加えることで,オブジェクトファイル(.o)が依存するソースファイル(.c, .cpp)とヘッダファイル(.h, .hpp)依存関係を自動的に定義することができる. あとは,実行ファイルとオブジェクトファイル間の依存関係を定義するだけでよい.

CXX := g++
SRCS := foo.cpp bar.cpp

$(foreach SRC,$(SRCS),$(eval $(subst \,,$(shell $(CXX) -MM $(SRC)))))

NTSC加重平均法の整数演算による高速化

はじめに

NTSC加重平均法という有名なグレースケール変換がある. このグレースケール変換は,ある画素の輝度値 $Y$ を,次の式(\ref{eq:ntscAverage})のように求めるものだ. ($r, g, b$ は赤,緑,青それぞれの値)

\begin{equation} Y = 0.298912 r + 0.586611 g + 0.114478 b \label{eq:ntscAverage} \end{equation}

単純な平均をとったグレースケール変換

\begin{equation} Y = \dfrac{r + g + b}{3} \end{equation}

と比較すると,NTSC加重平均法によるグレースケール変換は人間の赤,緑,青の感じ方を考慮したものになっており,自然なグレースケール化を行うことができる.

あるとき,このNTSC加重平均法によるグレースケール変換を,次のように整数のみで行っているプログラムを見ることがあった.

\begin{equation} Y = \dfrac{77 r + 150 g + 29 b}{256} \label{eq:ntscAverageOptimized08} \end{equation}

浮動小数点演算を整数演算に変換する発想に感動したので,この記事では,どのような過程で式(\ref{eq:ntscAverage})から式(\ref{eq:ntscAverageOptimized08})を得るのかについて解説する.

導出

記述の簡略化のため,

\begin{equation} W_r = 0.298912, \:\: W_g = 0.586611, \:\: W_b = 0.114478 \end{equation}

とおく. まず,理想的なNTSC加重平均法による輝度値を求める関数を以下の式(\ref{eq:ntscIdealFunction})ように定める.

\begin{equation} F(r, g, b) = W_r r + W_g g + W_b b \label{eq:ntscIdealFunction} \end{equation}

これを計算機で計算し,最終的に得られる輝度値(小数点を切り捨てたもの)を求める関数を以下の式(\ref{eq:ntscFlooredFunction})ように定める.

\begin{equation} F'(r, g, b) = \lfloor W_r r + W_g g + W_b b \rfloor \label{eq:ntscFlooredFunction} \end{equation}

ここで,関数 $F (r, g, b)$ を以下のように変形し, $f_N(r, g, b)$ を $F (r, g, b)$ の近似値を求める関数として定義する.

\begin{eqnarray} F(r, g, b) & = & F(r, g, b) \times \dfrac{2^N}{2^N} \nonumber \\ & = & \dfrac{2^N (W_r r + W_g g + W_b b)}{2^N} \nonumber \\ & \approx & \dfrac{\lfloor 2^N W_r + 0.5 \rfloor r + \lfloor 2^N W_g + 0.5 \rfloor g + \lfloor 2^N W_b + 0.5 \rfloor b}{2^N} \nonumber \\ & \equiv & f_N(r, g, b) \label{eq:ntscAverageApproxDerivation} \end{eqnarray}

$2^N$ を掛けているのは,2の累乗の計算はビット演算として表現可能なためである. ビット演算は掛け算や割り算と比べて,比較的高速に計算可能であり,容易に高速化が見込める.

さて, $N = 8$ のとき,すなわち, $f_8(r, g, b)$ は以下のようになる.

\begin{equation} f_8 (r, g, b) = \dfrac{77r + 150g + 29b}{256} \end{equation}

誤差についての考察

前述の導出より,式(\ref{eq:ntscAverageOptimized08})がNTSC加重平均法の近似式であることがわかった. しかし,近似計算に誤差はつきものである. ここでは,$N$ によって,特に $N = 8$ のとき,どの程度の誤差が生じるのかについて調べてみる.

まず,記述の簡略化のために,

\begin{equation} w_{rN} = 2^N \times W_r, \:\: w_{gN} = 2^N \times W_g, \:\: w_{bN} = 2^N \times W_b \end{equation}

\begin{equation} w'_{rN} = \lfloor w_{rN} + 0.5 \rfloor, \:\: w'_{gN} = \lfloor w_{gN} + 0.5 \rfloor, \:\: w'_{bN} = \lfloor w_{bN} + 0.5 \rfloor \end{equation}

とおく. 前述の導出式(\ref{eq:ntscAverageApproxDerivation})より,

\begin{equation} F (r, g, b) = f_N (r, g, b) + \dfrac{(w_{rN} - w'_{rN}) r + (w_{gN} - w'_{gN}) g + (w_{bN} - w'_{bN}) b}{2^N} \end{equation}

となる.右辺の第二項を理想の値との誤差を求める関数として,

\begin{equation} \epsilon_N (r, g, b) = \dfrac{(w_{rN} - w'_{rN}) r + (w_{gN} - w'_{gN}) g + (w_{bN} - w'_{bN}) b}{2^N} \end{equation}

とおく. ここに $N = 8$ を代入すると,

\begin{equation} \epsilon_8 (r, g, b) = \dfrac{-0.478528 r + 0.172416 g + 0.306368 b}{2^8} \end{equation}

を得る. $r, g, b$ の各係数の符号に注目すると,

\begin{equation} \epsilon_8 (255, 0, 0) \leq \epsilon_8 (r, g, b) \leq \epsilon_8 (0, 255, 255) \end{equation}

すなわち,

\begin{equation} -0.479 \leq \epsilon_8 (r, g, b) \leq 0.479 \end{equation}

を得る(小数第四位で四捨五入している). これは,やや誤差が大きいので,実際に計算機に計算させるときは, $N = 10$ として,

\begin{equation} \epsilon_{10} (0, 255, 0) = -0.0773 \leq \epsilon_{10} (r, g, b) \leq 0.0775 = \epsilon_{10} (255, 0, 255) \end{equation}

程度に誤差を抑えるのがよいだろう. もちろん, $N$ をより大きくすれば,誤差をさらに抑えることができる.そのときには,オーバーフローに気をつけなくてはならない.

ここで,r, g, bの各ビット深度が8bitであり,途中の計算結果を32bit整数として扱うと考える.

\begin{equation} W_r + W_g + W_b = 1 \end{equation}

かつ,

\begin{equation} 0 \leq r, g, b \leq 255 \end{equation}

であるので,

\begin{equation} 0 \leq W_r r + W_g g + W_b b < 2^8 \end{equation}

したがって,

\begin{equation} 0 \leq 2^{24} (W_r r + W_g g + W_b b) < 2^{32} \label{eq:ntscIdealShift} \end{equation}

となり, $N = 24$ がオーバーフローを起こさず,誤差を最も小さくできるビットシフト量である...と結論付けたいところであるが,実際には,$2^N W_r, \:\: 2^N W_g, \:\: 2^N W_b$ をそれぞれ四捨五入した係数 $w'_{rN}, \:\: w'_{gN}, \:\: w'_{bN}$ を用いるので,$N = 24$ のとき,式(\ref{eq:ntscIdealShift})が $2^{32}$ で抑えられるとは限らない. したがって, $N = 23$ が誤差を最も小さくできるビットシフト量であり,そのときの誤差 $\epsilon_{23} (r, g, b)$ は

\begin{equation} -2.061 \times 10^{-5} \leq \epsilon_{23} (r, g, b) \leq 2.025 \times 10^{-6} \end{equation}

となる.

まとめ

この記事では,NTSC加重平均法によるグレースケール変換を,整数演算のみで行う近似計算の導出を行い,その誤差について述べた.

具体的には,以下のような計算を行っているプログラム

y = (unsigned char) (0.298912 * r + 0.586611 * g + 0.114478 * b)

y = (unsigned char) ((306 * r + 601 * g + 117 + b) >> 10)

のように変更すると,元の計算とほとんど相違なく,高速に計算を行うことが可能となる.

グレースケール変換はOpenCVにも実装されているので,自分で実装する機会はほとんどないかもしれないが,整数 -> 整数という計算の途中に浮動小数点が介入する場合の高速化のアイデアとして,示唆に富むものがあるだろう.