$\phantom{}$$
\newcommand\nr{\notag\\}%タグなし改行用
\newcommand\ret{\notag\\&\qquad}%数式改行用
\newcommand\I{\mathrm{i}}
\newcommand\D{\mathrm{d}}
\newcommand\E{\mathrm{e}}
\newcommand\hc{\mathrm{h.c.}}
\newcommand\cc{\mathrm{c.c.}}
\newcommand\O[1]{\mathscr{O}\left(#1\right)}
\newcommand\abs[1]{{\left\rvert #1 \right\lvert}}
\newcommand\Res{\mathop{\mathrm{Res}}}
\newcommand\bra[1]{\mathinner{\langle{#1}|}}
\newcommand\ket[1]{\mathinner{|{#1}\rangle}}
\newcommand\braket[1]{\mathinner{\langle{#1}\rangle}}
\newcommand\Bra[1]{\left\langle#1\right|}
\newcommand\Ket[1]{\left|#1\right\rangle}
\newcommand\Braket[1]{\left\langle #1 \right\rangle}
\newcommand|{\middle|}%\Braket用
\DeclareMathOperator{\Log}{Log}
\DeclareMathOperator{\Arg}{Arg}
\DeclareMathOperator{\sgn}{sgn}
\DeclareMathOperator{\Tr}{Tr}
\newcommand\dashint{\mathchoice
{\int\kern-10pt-}
{\int\kern-8.5pt-}
{\int\kern-6.1pt-}
{\int\kern-4.58pt-}}
\newcommand\set[1]{\left\{#1\right\}}
\newcommand\for[1]{\quad\mathrm{for}\quad #1}
\newcommand\LHS{\mathrm{(LHS)}}
\newcommand\RHS{\mathrm{(RHS)}}
$
ACRiルームは，ACRiの活動の一貫として運営されている，100枚を超える FPGA ボードや FPGA StarterBOX を含むサーバ計算機をリモートからアクセスして利用できる FPGA 利用環境です．ACRi登録ページからアカウント申請をすれば使えるようになります．
最近，Intel(R) Agilex(TM) AGFB014R24B2E2V を搭載したFPGAボードDE10-Agilexが利用できるようになりました．使い方は，Agilexサーバーの利用方法を参考にしてもらうとして，この記事では，サンプルの中身を眺めながらDPC++で遊んでみることにします．題材として，量子コンピュータシミュレーションのQulacsのナイーブな移植を取り上げてみました．
サンプルを眺めてみる
ACRiルームの環境で利用できるoneAPI-sampleはoneAPI-samples 2022.3です．まずは，定番のTerasicのイントロダクションで取り上げられてる vector-add を見てみます．
一見長々としたコードですが，DPC++的な意味での主役は，VectorAdd 関数です．
void VectorAdd(queue &q, const IntVector &a_vector, const IntVector &b_vector,
IntVector &sum_parallel) {
// Create the range object for the vectors managed by the buffer.
range<1> num_items{a_vector.size()};
// Create buffers that hold the data shared between the host and the devices.
// The buffer destructor is responsible to copy the data back to host when it
// goes out of scope.
buffer a_buf(a_vector);
buffer b_buf(b_vector);
buffer sum_buf(sum_parallel.data(), num_items);
for (size_t i = 0; i < num_repetitions; i++ ) {
// Submit a command group to the queue by a lambda function that contains the
// data access permission and device computation (kernel).
q.submit([&](handler &h) {
// Create an accessor for each buffer with access permission: read, write or
// read/write. The accessor is a mean to access the memory in the buffer.
accessor a(a_buf, h, read_only);
accessor b(b_buf, h, read_only);
// The sum_accessor is used to store (with write permission) the sum data.
accessor sum(sum_buf, h, write_only, no_init);
// Use parallel_for to run vector addition in parallel on device. This
// executes the kernel.
// 1st parameter is the number of work items.
// 2nd parameter is the kernel, a lambda that specifies what to do per
// work item. The parameter of the lambda is the work item id.
// DPC++ supports unnamed lambda kernel by default.
h.parallel_for(num_items, [=](auto i) { sum[i] = a[i] + b[i]; });
});
};
// Wait until compute tasks on GPU done
q.wait();
}
forループは単に（？）num_repetitonsだけq.submitによるカーネル処理の実行を繰り返しているというだけで，処理の主役はh.parallel_forの部分です．ここで，num_items個の配列の要素をどかっと並列で足し算しています．
計算結果の確認のためにソフトウェアで計算されている部分を参照すると，DPC++を使わない場合の等価コードが確認できます．
for (size_t i = 0; i < sum_sequential.size(); i++)
sum_sequential.at(i) = a.at(i) + b.at(i);
の部分ですね．コードを読むと，num_items = a_vector.size() = sum_sequential.size()ということは見て取れます．
まねして遊んでみよう
題材選び
DPC++で遊ぶ雰囲気がわかったので，自分の書きたいプログラムを書いてみましょう．...と，ここで，なにやろうか？という話になるのですが，forループを単純に並列実行できる，しかもFPGAで実行できるということだけで遊べそうな，量子コンピュータのシミュレータをDPC++化してみましょう．
ここでは量子コンピュターの話に深く触れることはしませんが，量子コンピュータのシミュレータとして，n個の量子ビット列で表現されうる$2^n$個の複素数の確率振幅のベクトルに量子ゲート操作に相当するユニタリ行列をかけまくる，というモデルがあります．たとえば，Qulacsはそのようなシミュレータの高速なプログラムとして知られています．
$2^n$個の複素数の配列に対して行列をかけるだけなら，vector-addのサンプルと大差なさそうですね．というわけで，やってみましょう．
実装する処理の概要
簡単のため，量子コンピュータとして種々のアルゴリズムを動作させられる最低限の構成として，1量子ビットに対する操作と2量子に対する操作にだけ着目します．n量子ビットに対応する$n^2$個の複素数の配列に対して，前者は，操作を適用する量子ビットに対応する要素のペアに対して$2 \times 2$行列を掛けることに，後者は，操作を適用する量子ビットに対応する4つの要素セットに対して$ 4 \times 4$の行列を掛けることに相当します．
具体的には，1量子ビットの操作はupdate_ops_matrix_dense_single.cppのingle_qubit_dense_matrix_gate_single関数を2量子ビットの操作はupdate_ops_matrix_dense_doubleのdouble_qubit_dense_matrix_gate_nosimd関数をDPC++化することにします．
なお，コードを眺めるとわかるのですが，QulacsはSIMDやGPGPU向けにすごく最適化してあります．Qulacs: a fast and versatile quantum circuit simulator for research purposeはダテじゃないですね．
DPC++で実装してみる
題材が決まったので，実装してみましょう．
DPC++版 1量子ビット操作
というわけで，1量子ビット操作のコードを，parallel_forを使ってDPC++な並列化コードにします．オリジナルの関数の
for (state_index = 0; state_index < loop_dim; ++state_index) {
ITYPE basis_0 =
(state_index & mask_low) + ((state_index & mask_high) << 1);
ITYPE basis_1 = basis_0 + mask;
// fetch values
CTYPE cval_0 = state[basis_0];
CTYPE cval_1 = state[basis_1];
// set values
state[basis_0] = matrix[0] * cval_0 + matrix[1] * cval_1;
state[basis_1] = matrix[2] * cval_0 + matrix[3] * cval_1;
}
の中身を parallel_forに入れるラムダ関数化すれば良さそうですね．
というわけで，次のようにしてみました．型や変数名などはサンプルのvector_addに寄せています（寄せなければよかったとも思っています...）
//************************************
void single_qubit_gate(
queue &q,
int dim,
int target,
const std::vector<std::complex<double>> &U_vector,
const std::vector<std::complex<double>> &psi_vector,
std::vector<std::complex<double>> &psi_out_vector){
const unsigned int target_mask = 1 << target;
const unsigned int loop_dim = dim / 2;
range<1> num_items{psi_vector.size()};
range<1> num_loop_dim{loop_dim};
buffer U_buf(U_vector);
buffer psi_buf(psi_vector);
buffer psi_out_buf(psi_out_vector.data(), num_items);
q.submit([&](handler &h) {
accessor U(U_buf, h, read_only);
accessor psi(psi_buf, h, read_only);
accessor psi_out(psi_out_buf, h, write_only, no_init);
h.parallel_for(num_loop_dim, [=](auto i) {
int temp_basis = (i >> target) << (target + 1);
int basis_0 = temp_basis + i % target_mask;
int basis_1 = basis_0 ^ target_mask;
std::complex<double> cval_0 = psi[basis_0];
std::complex<double> cval_1 = psi[basis_1];
psi_out[basis_0] = U[0] * cval_0 + U[1] * cval_1;
psi_out[basis_1] = U[2] * cval_0 + U[3] * cval_1;
});
});
q.wait();
}
DPC++のサンプルvector_addでは，足し算の対象と成るベクトルの要素数をparallel_forのサイズに与えています．しかしながら，移植元のコードでは一度に2要素ずつの計算を行うために，parallel_forに与えるサイズloop_dimが対象となる配列の半分であることに注意する必要があります．まあ，forループの繰り返し回数をそのまま指定すればいい，というわけです．
2量子ビット操作の移植
2量子ビットの移植も同様です．Qulacsで，
for (state_index = 0; state_index < loop_dim; ++state_index) {
// create index
ITYPE basis_0 = (state_index & low_mask) +
((state_index & mid_mask) << 1) +
((state_index & high_mask) << 2);
// gather index
ITYPE basis_1 = basis_0 + target_mask1;
ITYPE basis_2 = basis_0 + target_mask2;
ITYPE basis_3 = basis_1 + target_mask2;
// fetch values
CTYPE cval_0 = state[basis_0];
CTYPE cval_1 = state[basis_1];
CTYPE cval_2 = state[basis_2];
CTYPE cval_3 = state[basis_3];
// set values
state[basis_0] = matrix[0] * cval_0 + matrix[1] * cval_1 +
matrix[2] * cval_2 + matrix[3] * cval_3;
state[basis_1] = matrix[4] * cval_0 + matrix[5] * cval_1 +
matrix[6] * cval_2 + matrix[7] * cval_3;
state[basis_2] = matrix[8] * cval_0 + matrix[9] * cval_1 +
matrix[10] * cval_2 + matrix[11] * cval_3;
state[basis_3] = matrix[12] * cval_0 + matrix[13] * cval_1 +
matrix[14] * cval_2 + matrix[15] * cval_3;
}
とループを回している部分を，DPC++のparallel_forを使って並列化します．こんな感じにしてみました．
void double_qubit_dense_matrix_gate_single(
queue &q,
int dim,
int target_qubit_index1,
int target_qubit_index2,
const std::vector<std::complex<double>> U_vector,
const std::vector<std::complex<double>> &psi_vector,
std::vector<std::complex<double>> &psi_out_vector
){
const unsigned int min_qubit_index = std::min(target_qubit_index1, target_qubit_index2);
const unsigned int max_qubit_index = std::max(target_qubit_index1, target_qubit_index2);
const unsigned int min_qubit_mask = 1ULL << min_qubit_index;
const unsigned int max_qubit_mask = 1ULL << (max_qubit_index - 1);
const unsigned int low_mask = min_qubit_mask - 1;
const unsigned int mid_mask = (max_qubit_mask - 1) ^ low_mask;
const unsigned int high_mask = ~(max_qubit_mask - 1);
const unsigned int target_mask1 = 1ULL << target_qubit_index1;
const unsigned int target_mask2 = 1ULL << target_qubit_index2;
// loop variables
const unsigned int loop_dim = dim / 4;
range<1> num_items{psi_vector.size()};
range<1> num_loop_dim{loop_dim};
buffer U_buf(U_vector);
buffer psi_buf(psi_vector);
buffer psi_out_buf(psi_out_vector.data(), num_items);
q.submit([&](handler &h) {
accessor U(U_buf, h, read_only);
accessor psi(psi_buf, h, read_only);
accessor psi_out(psi_out_buf, h, write_only, no_init);
h.parallel_for(num_loop_dim, [=](auto i) {
// create index
unsigned int basis_0 = (i & low_mask) + ((i & mid_mask) << 1) + ((i & high_mask) << 2);
// gather index
unsigned int basis_1 = basis_0 + target_mask1;
unsigned int basis_2 = basis_0 + target_mask2;
unsigned int basis_3 = basis_1 + target_mask2;
// fetch values
std::complex<double> cval_0 = psi[basis_0];
std::complex<double> cval_1 = psi[basis_1];
std::complex<double> cval_2 = psi[basis_2];
std::complex<double> cval_3 = psi[basis_3];
// set values
psi_out[basis_0] = U[0] * cval_0 + U[1] * cval_1 + U[2] * cval_2 + U[3] * cval_3;
psi_out[basis_1] = U[4] * cval_0 + U[5] * cval_1 + U[6] * cval_2 + U[7] * cval_3;
psi_out[basis_2] = U[8] * cval_0 + U[9] * cval_1 + U[10] * cval_2 + U[11] * cval_3;
psi_out[basis_3] = U[12] * cval_0 + U[13] * cval_1 + U[14] * cval_2 + U[15] * cval_3;
});
});
q.wait();
}
こちらも，ループ部分を素直にparallel_forに置き換えただけ，です．
動作確認してみる
1量子ビットと2量子ビットに相当する関数が用意できたので，動作確認をしてみましょう．代表的な1量子ビットと2量子ビットの演算にアダマール演算とCNOT演算があります．それぞれ，$2 \times 2$あるいは$4 \times 4$に相当する行列として，次のように定義できます．
std::vector<std::complex<double>> H;
H.push_back(std::complex<double>(1.0/sqrt(2), 0.0f));
H.push_back(std::complex<double>(1.0/sqrt(2), 0.0f));
H.push_back(std::complex<double>(1.0/sqrt(2), 0.0f));
H.push_back(std::complex<double>(-1.0/sqrt(2), 0.0f));
std::vector<std::complex<double>> CNOT;
CNOT.push_back(std::complex<double>(1.0, 0.0f));
CNOT.push_back(std::complex<double>(0.0, 0.0f));
CNOT.push_back(std::complex<double>(0.0, 0.0f));
CNOT.push_back(std::complex<double>(0.0, 0.0f));
CNOT.push_back(std::complex<double>(0.0, 0.0f));
CNOT.push_back(std::complex<double>(1.0, 0.0f));
CNOT.push_back(std::complex<double>(0.0, 0.0f));
CNOT.push_back(std::complex<double>(0.0, 0.0f));
CNOT.push_back(std::complex<double>(0.0, 0.0f));
CNOT.push_back(std::complex<double>(0.0, 0.0f));
CNOT.push_back(std::complex<double>(0.0, 0.0f));
CNOT.push_back(std::complex<double>(1.0, 0.0f));
CNOT.push_back(std::complex<double>(0.0, 0.0f));
CNOT.push_back(std::complex<double>(0.0, 0.0f));
CNOT.push_back(std::complex<double>(1.0, 0.0f));
CNOT.push_back(std::complex<double>(0.0, 0.0f));
このように量子演算に相当する行列（C++的にはstd::vector）を用意しておいて，
queue q(d_selector, exception_handler);
// Print out the device information used for the kernel code.
std::cout << "Running on device: "
<< q.get_device().get_info<info::device::name>() << std::endl;
// Vector addition in DPC++
single_qubit_gate(q, dim, target_index, H, psi, psi_out);
print_state(psi);
std::cout << "H" << std::endl;
print_state(psi_out);
double_qubit_dense_matrix_gate_single(q, dim, cnot_control_index, cnot_target_index, CNOT, psi_out, psi);
std::cout << "CNOT" << std::endl;
print_state(psi);
などとすることで，実装したDPC++版のコードを使った量子コンピュータシミュレーションを試すことができます．
ここでは，n個の量子ビット(複素数の配列サイズは$dim = 2^n$)に対して，target_indexで与えた量子ビットにアダマール演算を適用し，そのあとで，cnot_control_indexとcnot_target_indexで指定したCNOT演算を適用しています．
実行してみた
わざわざFPGAを使わなくてもソフトウェア上で動作確認ができるのもDPC++の嬉しさです．
miyo@iserv1:~/oneAPI-samples-2022.3.0/DirectProgramming/DPC++/DenseLinearAlgebra/vector-add$ ./vector-add-buffers 3 0 1 0
qubit=3, dim=8
H: target=0
CNOT: control=1, target=0
Running on device: 11th Gen Intel(R) Core(TM) i7-11700K @ 3.60GHz
[0]: (1,0)
[1]: (0,0)
[2]: (0,0)
[3]: (0,0)
[4]: (0,0)
[5]: (0,0)
[6]: (0,0)
[7]: (0,0)
H
[0]: (0.707107,0)
[1]: (0.707107,0)
[2]: (0,0)
[3]: (0,0)
[4]: (0,0)
[5]: (0,0)
[6]: (0,0)
[7]: (0,0)
CNOT
[0]: (0.707107,0)
[1]: (0,0)
[2]: (0,0)
[3]: (0.707107,0)
[4]: (0,0)
[5]: (0,0)
[6]: (0,0)
[7]: (0,0)
miyo@iserv1:~/oneAPI-samples-2022.3.0/DirectProgramming/DPC++/DenseLinearAlgebra/vector-add$
といった具合に，たとえば3量子ビットがある場合に，初期状態 $\ket{000}$ (なので，[0]番目の確率振幅だけが1になっている）にアダマール演算を適用して重ね合わせ状態を作って，CNOTで伝播させた様子が見て取れます．
FPGA向けにビルドして実行する
もちろんFPGA向けに合成して実行することもできます．vector-addサンプルからスタートしたので，Makefileをそのまま利用することができます．コンパイルのコマンドは具体的には次の通りです．
CXX := dpcpp
CXXFLAGS = -O2 -g -std=c++17
SRC := src/vector-add-buffers.cpp
hw: vector-add-buffers.fpga
a.o: $(SRC)
$(CXX) $(CXXFLAGS) -fintelfpga -c $^ -o $@ -DFPGA=1
vector-add-buffers.fpga: a.o
$(CXX) $(CXXFLAGS) -fintelfpga $^ -o $@ -Xshardware -Xsboard=de10_agilex:B2E2_8GBx4
ビルドした様子は次の通り.
dpcpp -O2 -g -std=c++17 -fintelfpga -c src/vector-add-buffers.cpp -o a.o -DFPGA=1
dpcpp -O2 -g -std=c++17 -fintelfpga a.o -o vector-add-buffers.fpga -Xshardware -Xsboard=de10_agilex:B2E2_8GBx4
aoc: Compiling for FPGA. This process may take several hours to complete. Prior to performing this compile, be sure to check the reports to ensure the design will meet your performance targets. If the reports indicate performance targets are not being met, code edits may be required. Please refer to the oneAPI FPGA Optimization Guide for information on performance tuning applications for FPGAs.
'quartus_compile_report.log' and '/tmp/a-5ecfd7-805af9/quartus_compile_report.log' are identical (not copied) at /opt/intel/oneapi/compiler/2022.0.2/linux/lib/oclfpga/share/lib/perl/acl/aoc.pl line 1578.
ビルド結果を見てみると
itter Status : Successful - Sat Dec 24 20:13:37 2022
Quartus Prime Version : 21.2.0 Build 72 06/14/2021 SC Pro Edition
Revision Name : top
Top-level Entity Name : top
Family : Agilex
Device : AGFB014R24B2E2V
Timing Models : Preliminary
Power Models : Preliminary
Device Status : Preliminary
Logic utilization (in ALMs) : 179,945 / 487,200 ( 37 % )
Total dedicated logic registers : 504415
Total pins : 567 / 924 ( 61 % )
Total block memory bits : 18,069,984 / 145,612,800 ( 12 % )
Total RAM Blocks : 1,366 / 7,110 ( 19 % )
Total DSP Blocks : 240 / 4,510 ( 5 % )
Total eSRAMs : 0 / 2 ( 0 % )
Total HSSI P-Tiles : 1 / 1 ( 100 % )
Total HSSI E-Tile Channels : 0 / 24 ( 0 % )
Total HSSI HPS : 0 / 1 ( 0 % )
Total HSSI EHIPs : 0 / 4 ( 0 % )
Total PLLs : 16 / 24 ( 67 % )
でした．実行結果はこんな感じ．ソフトウェアでの実行と同じですね．
miyo@iserv1:~/oneAPI-samples-2022.3.0/DirectProgramming/DPC++/DenseLinearAlgebra/vector-add$ ./vector-add-buffers.fpga 3 0 1 0
qubit=3, dim=8
H: target=0
CNOT: control=1, target=0
Running on device: B2E2_8GBx4 : Agilex Reference Platform (aclde10_agilex0)
[0]: (1,0)
[1]: (0,0)
[2]: (0,0)
[3]: (0,0)
[4]: (0,0)
[5]: (0,0)
[6]: (0,0)
[7]: (0,0)
H
[0]: (0.707107,0)
[1]: (0.707107,0)
[2]: (0,0)
[3]: (0,0)
[4]: (0,0)
[5]: (0,0)
[6]: (0,0)
[7]: (0,0)
CNOT
[0]: (0.707107,0)
[1]: (0,0)
[2]: (0,0)
[3]: (0.707107,0)
[4]: (0,0)
[5]: (0,0)
[6]: (0,0)
[7]: (0,0)
miyo@iserv1:~/oneAPI-samples-2022.3.0/DirectProgramming/DPC++/DenseLinearAlgebra/vector-add$
実用的に（？）使うには
さすがに，シミュレーションしたい量子プログラムを適用するためにC++ソースコードを変更してコンパイル（内部ではQuartusによるビルドが走る！）を実行するのはバカバカしいですね．というわけで，何かしらのフロントエンドとの連携は必要そうです．とはいえ，Qulacs自体が十分高速なので，DPC++/FPGAを使った価値を見出すための高速化に先に取り組む方が良さそうです．
というわけで
この記事では，vector-addを紹介して，とりあえず手軽に既存のC++のループをDPC++のparallel_forにおきかえてみるという遊びができますよ，という話をしてみました．oneAPI-samplesに他にもたくさんありますので，自分でコード書く前にもっとサンプルみたいという場合にも，是非いろいろ試してみてください．
DPC++楽しそう，ACRiルームならFPGA使うにせよ使わないにせよ構築済みの環境で手軽に試せそう，と思っていただければ嬉しいです．...あれ，こんな終わり方でいいんだっけ？