去年のアドベントカレンダーで量子スピン模型のハミルトニアンを数値対角化しようというタイトルで記事を書いたのですが、やっぱり今時はGPUだよね!ということで、去年のコードをGPUで走らせてみたいと思います。
ところで最近monolishというライブラリがC++から疎行列計算が簡単に扱えると聞いたので、このライブラリを使ってみたいと思います。
なので、今年のコードはC++です。
ちなみに、この記事は最初は副題の通りmonolishの使い方ガイドとして書いてたので、そんな感じになってます。
monolishってなぁに
すべてのハードウェア、型、行列構造への依存を吸収することを目的として開発された、数値計算ライブラリ"monolish"を公開しました。https://t.co/g5ZLJGAlGy
— 科学計算総合研究所 (RICOS Co.Ltd.) (@RICOS_ltd) April 5, 2021
PythonやJuliaのような簡単なインタフェースで、MKL、 cuBLAS、cuSPARSE、cuSOLVER、BLAS、LAPACKのAPIを統一した関数が利用できます。
monolishの最も気持ち悪いコードはこれです
— 電子計算機の沼 (@Hishinuma_t) April 5, 2021
俺がクソみたいな全ての沼を飲み干してやると言うつもりで作りましたhttps://t.co/9FKPJGUymF
私にとってはC++でCPU版とGPU版のコードが共通化できるライブラリという認識です。
2021年12月8日現在、ソースコードはこちらでドキュメントはこちらです。
実行環境
Ubuntu 18.04.01 LTS
CPU: Intel(R) Core(TM) i5-8400
GPU: GeForce RTX 2070
CUDA Version: 11.2
monolish: version 0.14.0
monolishをインストールする
この記事では提供されているdockerを使います。
docker慣れてないですけど頑張りました。
dockerのinstallにはこちらの記事を参考にしました。
monolishのgithubページからソースコードをcloneしてきます。
git clone https://github.com/ricosjp/monolish.git
monolishの中に入ります
cd monolish
そこでmake in_oss_cpu
(GPU版の場合はmake in_oss_gpu
)と打ち込みます。
make in_oss_cpu
(または make in_oss_gpu)
すると勝手にdockerの中に入ります。
gpu版の場合、dockerに入ろうとするタイミングで
docker: Error response from daemon: linux runtime spec devices: could not select device driver "" with capabilities: [[gpu]].
のエラーが出るかもしれません。どうやらdockerでgpuを動かす時のエラーのようです。私の場合はgoogle検索で出てきたサイト(https://github.com/NVIDIA/nvidia-docker/issues/1034 とか)を片っ端から試して全部ダメでふて寝して起きてもう一度やろうとしたらdockerコマンドが壊れてたので再インストールしたら出来るようになりました。何が原因かわかりません。
dockerに入ったらそこでmake isntall_cpu
またはmake install_gpu
をします。
make install_cpu
(または make install_gpu)
これでmonolishを使える環境が整います。
monolishのサンプルコード
examples内(もしくはgithubのページ)にサンプルコードがあります。
やってみる
量子スピン模型のハミルトニアンを数値対角化しようのコードをmonolishを使ったC++にしたものがこちらになります。
#include<cstdio>
#include<monolish_eigen.hpp>
#include<monolish_equation.hpp>
#include<vector>
#include<complex>
#include<cmath>
#include<chrono>
const std::complex<double> Jx(1.0, 0.0);
const std::complex<double> Jy(1.0, 0.0);
const std::complex<double> Jz(1.0, 0.0);
typedef std::vector<std::vector<std::complex<double>>> OPERATOR;
void set_Hamiltonian_tmp(const OPERATOR O1, const int index1, const OPERATOR O2, const int index2, const int x, const int y, const int itr, const std::complex<double>& val, std::vector<int>& rows, std::vector<int>& cols, std::vector<double>& vals){
if(itr>=0){
if(itr==index1){
set_Hamiltonian_tmp(O1, index1, O2, index2, 2*x, 2*y, itr-1, val*O1[0][0], rows, cols, vals);
set_Hamiltonian_tmp(O1, index1, O2, index2, 2*x+1, 2*y, itr-1, val*O1[1][0], rows, cols, vals);
set_Hamiltonian_tmp(O1, index1, O2, index2, 2*x, 2*y+1, itr-1, val*O1[0][1], rows, cols, vals);
set_Hamiltonian_tmp(O1, index1, O2, index2, 2*x+1, 2*y+1, itr-1, val*O1[1][1], rows, cols, vals);
}else if(itr==index2){
set_Hamiltonian_tmp(O1, index1, O2, index2, 2*x, 2*y, itr-1, val*O2[0][0], rows, cols, vals);
set_Hamiltonian_tmp(O1, index1, O2, index2, 2*x+1, 2*y, itr-1, val*O2[1][0], rows, cols, vals);
set_Hamiltonian_tmp(O1, index1, O2, index2, 2*x, 2*y+1, itr-1, val*O2[0][1], rows, cols, vals);
set_Hamiltonian_tmp(O1, index1, O2, index2, 2*x+1, 2*y+1, itr-1, val*O2[1][1], rows, cols, vals);
}else{
set_Hamiltonian_tmp(O1, index1, O2, index2, 2*x, 2*y, itr-1, val, rows, cols, vals);
set_Hamiltonian_tmp(O1, index1, O2, index2, 2*x+1, 2*y+1, itr-1, val, rows, cols, vals);
}
}else{
rows.push_back(x);
cols.push_back(y);
vals.push_back(val.real());
}
return;
}
void set_Hamiltonian(const OPERATOR O1, int index1, const OPERATOR O2, int index2, const std::complex<double>& J, std::vector<int>& rows, std::vector<int>& cols, std::vector<double>& vals, const int N){
set_Hamiltonian_tmp(O1, index1, O2, index2, 0, 0, N-1, J, rows, cols, vals);
return;
}
std::vector<std::complex<double>> set_Hamiltonian_diag(const OPERATOR O1, int index1, const OPERATOR O2, int index2, const std::complex<double>& J, const int N){
std::vector<std::complex<double>> val_array(1, J);
for(int i=N-1; i>=0; --i){
val_array.resize(val_array.size()*2);
if(i==index1){
for(long long j=val_array.size()/2-1; j>=0; --j){
val_array[2*j+1] = val_array[j]*O1[1][1];
val_array[2*j] = val_array[j]*O1[0][0];
}
}else if(i==index2){
for(long long j=val_array.size()/2-1; j>=0; --j){
val_array[2*j+1] = val_array[j]*O2[1][1];
val_array[2*j] = val_array[j]*O2[0][0];
}
}else{
for(long long j=val_array.size()/2-1; j>=0; --j){
val_array[2*j+1] = val_array[j];
val_array[2*j] = val_array[j];
}
}
}
return val_array;
}
void set_Hamiltonian_offdiag(const OPERATOR& O1_1, const OPERATOR& O2_1, const int index1, const OPERATOR& O1_2, const OPERATOR& O2_2, const int index2, const std::complex<double>& J1, const std::complex<double>& J2, const int N, std::vector<int>& rows, std::vector<int>& cols, std::vector<double>& vals){
std::vector<int> x_array(1, 0);
std::vector<int> y_array(1, 0);
std::vector<std::complex<double>> val1_array(1, J1);
std::vector<std::complex<double>> val2_array(1, J2);
for(int i=N-1; i>=0; --i){
x_array.resize(x_array.size()*2);
y_array.resize(y_array.size()*2);
val1_array.resize(val1_array.size()*2);
val2_array.resize(val2_array.size()*2);
if(i==index1){
for(long long j=x_array.size()/2-1; j>=0; --j){
x_array[2*j+1] = 2*x_array[j]+1;
y_array[2*j+1] = 2*y_array[j];
val1_array[2*j+1] = val1_array[j]*O1_1[1][0];
val2_array[2*j+1] = val2_array[j]*O2_1[1][0];
x_array[2*j] = 2*x_array[j];
y_array[2*j] = 2*y_array[j]+1;
val1_array[2*j] = val1_array[j]*O1_1[0][1];
val2_array[2*j] = val2_array[j]*O2_1[0][1];
}
}else if(i==index2){
for(long long j=x_array.size()/2-1; j>=0; --j){
x_array[2*j+1] = 2*x_array[j]+1;
y_array[2*j+1] = 2*y_array[j];
val1_array[2*j+1] = val1_array[j]*O1_2[1][0];
val2_array[2*j+1] = val2_array[j]*O2_2[1][0];
x_array[2*j] = 2*x_array[j];
y_array[2*j] = 2*y_array[j]+1;
val1_array[2*j] = val1_array[j]*O1_2[0][1];
val2_array[2*j] = val2_array[j]*O2_2[0][1];
}
}else{
for(long long j=x_array.size()/2-1; j>=0; --j){
x_array[2*j+1] = 2*x_array[j]+1;
y_array[2*j+1] = 2*y_array[j]+1;
val1_array[2*j+1] = val1_array[j];
val2_array[2*j+1] = val2_array[j];
x_array[2*j] = 2*x_array[j];
y_array[2*j] = 2*y_array[j];
val1_array[2*j] = val1_array[j];
val2_array[2*j] = val2_array[j];
}
}
}
rows.insert(rows.end(), x_array.begin(), x_array.end());
cols.insert(cols.end(), y_array.begin(), y_array.end());
size_t vals_size = vals.size();
vals.resize(vals.size()+val1_array.size());
for(size_t i=0; i<val1_array.size(); ++i){
vals[vals_size+i] = (val1_array[i]+val2_array[i]).real();
}
return;
}
void initialize(const int N, std::vector<int>& rows, std::vector<int>& cols, std::vector<double>& vals,
const OPERATOR& Sx, const OPERATOR& Sy, const OPERATOR& Sz
){
for(int i=0; i<N-1; ++i){
set_Hamiltonian(Sx, i, Sx, i+1, Jx, rows, cols, vals, N);
set_Hamiltonian(Sy, i, Sy, i+1, Jy, rows, cols, vals, N);
set_Hamiltonian(Sz, i, Sz, i+1, Jz, rows, cols, vals, N);
}
set_Hamiltonian(Sx, N-1, Sx, 0, Jx, rows, cols, vals, N);
set_Hamiltonian(Sy, N-1, Sy, 0, Jy, rows, cols, vals, N);
set_Hamiltonian(Sz, N-1, Sz, 0, Jz, rows, cols, vals, N);
return;
}
void initialize_diag(const size_t N, std::vector<int>& rows, std::vector<int>& cols, std::vector<double>& vals, const OPERATOR& Sx, const OPERATOR& Sy, const OPERATOR& Sz){
size_t M = pow(2, N);
std::vector<std::complex<double>> diag(M);
for(size_t i=0; i<N-1; ++i){
std::vector<std::complex<double>> diag_tmp = set_Hamiltonian_diag(Sz, i, Sz, i+1, Jz, N);
for(size_t j=0; j<diag.size(); ++j){
diag[j] += diag_tmp[j];
}
}
std::vector<std::complex<double>> diag_tmp = set_Hamiltonian_diag(Sz, N-1, Sz, 0, Jz, N);
for(size_t j=0; j<diag.size(); ++j){
diag[j] += diag_tmp[j];
}
for(size_t j=0; j<diag.size(); ++j){
rows.push_back(j);
cols.push_back(j);
vals.push_back(diag[j].real());
}
return;
}
void initialize_offdiag(const size_t N, std::vector<int>& rows, std::vector<int>& cols, std::vector<double>& vals, const OPERATOR& Sx, const OPERATOR& Sy, const OPERATOR& Sz){
for(size_t i=0; i<N-1; ++i){
set_Hamiltonian_offdiag(Sx, Sy, i, Sx, Sy, i+1, Jx, Jy, N, rows, cols, vals);
}
set_Hamiltonian_offdiag(Sx, Sy, N-1, Sx, Sy, 0, Jx, Jy, N, rows, cols, vals);
return;
}
void convert(const std::vector<int>& rows, const std::vector<int>& cols, const std::vector<double>& vals, std::vector<int>& rowptr, std::vector<int>& colvec, std::vector<double>& valvec){
std::vector<std::vector<int>> cols_tmp(rowptr.size()-1);
std::vector<std::vector<double>> vals_tmp(rowptr.size()-1);
if(rows.size() != cols.size()){
throw std::runtime_error("rows.size != cols.size");
}
if(rows.size() != vals.size()){
throw std::runtime_error("rows.size != vals.size");
}
for(size_t i=0; i<rows.size(); ++i){
cols_tmp[rows[i]].push_back(cols[i]);
vals_tmp[rows[i]].push_back(vals[i]);
}
int begin = 0;
for(size_t i=0; i<cols_tmp.size(); ++i){
rowptr[i] = begin;
colvec.insert(colvec.end(), cols_tmp[i].begin(), cols_tmp[i].end());
valvec.insert(valvec.end(), vals_tmp[i].begin(), vals_tmp[i].end());
begin += cols_tmp[i].size();
}
rowptr[cols_tmp.size()] = begin;
return;
}
int main(int argc, char** argv){
if(argc<2){
fprintf(stderr, "usage : %s N\n", argv[0]);
exit(0);
}
std::chrono::system_clock::time_point start, end;
start = std::chrono::system_clock::now();
size_t N = std::atoi(argv[1]);
size_t M = pow(2, N);
size_t K = 6;
fprintf(stderr, "N=%zu\n", N);
OPERATOR Sx(2, std::vector<std::complex<double>>(2));
OPERATOR Sy(2, std::vector<std::complex<double>>(2));
OPERATOR Sz(2, std::vector<std::complex<double>>(2));
Sx[0][1] = 1; Sx[1][0] = 1;
Sy[0][1] = std::complex<double>(0, -1); Sy[1][0] = std::complex<double>(0, 1);
Sz[0][0] = 1; Sz[1][1] = -1;
std::vector<int> rows, cols;
std::vector<double> vals;
initialize_offdiag(N, rows, cols, vals, Sx, Sy, Sz);
initialize_diag(N, rows, cols, vals, Sx, Sy, Sz);
std::vector<int> rowptr(M+1), colvec;
std::vector<double> valvec;
convert(rows, cols, vals, rowptr, colvec, valvec);
monolish::matrix::CRS<double> MAT(M, M, rowptr, colvec, valvec);
monolish::vector<double> lambda(K);
monolish::matrix::Dense<double> x(K, MAT.get_col());
monolish::eigen::LOBPCG<monolish::matrix::CRS<double>, double> solver;
//monolish::equation::Jacobi<monolish::matrix::CRS<double>, double> precond;
//solver.set_create_precond(precond);
//solver.set_apply_precond(precond);
double tol_res = 1.0e-8;
solver.set_tol(tol_res);
solver.set_lib(0);
solver.set_miniter(0);
int maxiter = N*10;
solver.set_maxiter(maxiter);
end = std::chrono::system_clock::now();
fprintf(stderr, "initialize time = %lf\n", ((double)std::chrono::duration_cast<std::chrono::milliseconds>(end-start).count())/1000);
fprintf(stderr, "eig start\n");
start = std::chrono::system_clock::now();
monolish::util::send(MAT);
monolish::util::solver_check(solver.solve(MAT, lambda, x));
monolish::util::recv(lambda, x);
end = std::chrono::system_clock::now();
fprintf(stderr, "eig time = %lf\n", ((double)std::chrono::duration_cast<std::chrono::milliseconds>(end-start).count())/1000);
fprintf(stderr, "eig end\n");
fprintf(stderr, "%.16lf\n", lambda[0]/N);
return 0;
}
量子スピン模型のハミルトニアンを数値対角化しようの時とは違い、CRS形式の疎行列に渡す引数はCRS形式でないとダメなので、main
関数内で行列の情報rows, cols, vals(疎行列Matに対してMat[rows[i]][cols[i]] = vals[i]となっている)を生成した後、converter
関数でrows, cols, vals を rowptr, colvec, valves に変換しています。
CRS形式の引数としてはwikipediaのページなどを参考にしました。
上記ソースコード内の monolish と書いてるところが今回のmonolish関係のコードです。サンプルコード通りにprecondを行うと答えが合わなかったのでコメントアウトしてあります。
send
/recv
はGPU版で必要(GPUにデータを送る/GPUからデータを返す)で、CPU版では必要はないです。書いてても動きます。
g++ spin.cpp -o spin -O3 -I /opt/monolish/include -L /opt/monolish/lib -lmonolish_cpu
とするとコンパイルできます。(GPU版の場合は-lmonolish_gpu
)
リンクを忘れずに!!!
計算結果比較
N=16で以前のpythonの結果が-1.7855740901541892、今回の結果が-1.7855740901541948なので、だいたい合ってそうです
計算時間比較
計算時間の比較が次の図です。GPUは転送時間を含みます。
ちゃんとGPUが使えてて、N=22ではGPUを使って3倍高速化できてそうです。
N=16などの小さいところでは、GPUとの通信のせいでCPUよりも計算が遅くなっていそうです。
N=24はGPUのメモリが足りませんでした。。。
メモリたくさんのGPU欲しい。。。
monolishのインストールに半年で4回計24時間ぐらい失敗し続けたので、この記事が多少の助けになると幸いです。