スピン模型をGPUを使って解く　（副題：monolishを使って固有値を計算してみる）



去年のアドベントカレンダーで[量子スピン模型のハミルトニアンを数値対角化しよう](https://qiita.com/fockl/items/52593901b9a5f78cda92)というタイトルで記事を書いたのですが、やっぱり今時はGPUだよね！ということで、去年のコードをGPUで走らせてみたいと思います。

ところで最近monolishというライブラリがC++から疎行列計算が簡単に扱えると聞いたので、このライブラリを使ってみたいと思います。
なので、今年のコードはC++です。

ちなみに、この記事は最初は副題の通りmonolishの使い方ガイドとして書いてたので、そんな感じになってます。

# monolishってなぁに

<blockquote class="twitter-tweet" data-partner="tweetdeck"><p lang="ja" dir="ltr">すべてのハードウェア、型、行列構造への依存を吸収することを目的として開発された、数値計算ライブラリ&quot;monolish&quot;を公開しました。<a href="https://t.co/g5ZLJGAlGy">https://t.co/g5ZLJGAlGy</a><br><br>PythonやJuliaのような簡単なインタフェースで、MKL、 cuBLAS、cuSPARSE、cuSOLVER、BLAS、LAPACKのAPIを統一した関数が利用できます。</p>&mdash; 科学計算総合研究所 (RICOS Co.Ltd.) (@RICOS_ltd) <a href="https://twitter.com/RICOS_ltd/status/1379038168834400262?ref_src=twsrc%5Etfw">April 5, 2021</a></blockquote>
<script async src="https://platform.twitter.com/widgets.js" charset="utf-8"></script>

<blockquote class="twitter-tweet" data-partner="tweetdeck"><p lang="ja" dir="ltr">monolishの最も気持ち悪いコードはこれです<br><br>俺がクソみたいな全ての沼を飲み干してやると言うつもりで作りました<a href="https://t.co/9FKPJGUymF">https://t.co/9FKPJGUymF</a></p>&mdash; 電子計算機の沼 (@Hishinuma_t) <a href="https://twitter.com/Hishinuma_t/status/1379039625717587971?ref_src=twsrc%5Etfw">April 5, 2021</a></blockquote>
<script async src="https://platform.twitter.com/widgets.js" charset="utf-8"></script>

私にとってはC++でCPU版とGPU版のコードが共通化できるライブラリという認識です。
2021年12月8日現在、ソースコードは[こちら](https://github.com/ricosjp/monolish)でドキュメントは[こちら](https://ricosjp.github.io/monolish/0.15.0/)です。


# 実行環境

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には[こちらの記事](https://qiita.com/iganari/items/fe4889943f22fd63692a)を参考にしました。

[monolishのgithubページ](https://github.com/ricosjp/monolish/)からソースコードをcloneしてきます。

```bash
git clone https://github.com/ricosjp/monolish.git
```

monolishの中に入ります

```bash
cd monolish
```

そこで`make in_oss_cpu` (GPU版の場合は`make in_oss_gpu`)と打ち込みます。

```bash
make in_oss_cpu
(または make in_oss_gpu)
```

すると勝手にdockerの中に入ります。
gpu版の場合、dockerに入ろうとするタイミングで

```bash
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`をします。

```bash
make install_cpu
(または make install_gpu)
```

これでmonolishを使える環境が整います。

# monolishのサンプルコード

examples内(もしくは[githubのページ](https://github.com/ricosjp/monolish/tree/master/examples))にサンプルコードがあります。

# やってみる

[量子スピン模型のハミルトニアンを数値対角化しよう](https://qiita.com/fockl/items/52593901b9a5f78cda92)のコードをmonolishを使ったC++にしたものがこちらになります。

```c++:spin.cpp
#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;
}
```

[量子スピン模型のハミルトニアンを数値対角化しよう](https://qiita.com/fockl/items/52593901b9a5f78cda92)の時とは違い、CRS形式の疎行列に渡す引数はCRS形式でないとダメなので、`main`関数内で行列の情報rows, cols, vals（疎行列Matに対してMat[rows[i]][cols[i]] = vals[i]となっている）を生成した後、`converter`関数でrows, cols, vals を rowptr, colvec, valves に変換しています。
CRS形式の引数としては[wikipediaのページ](https://ja.wikipedia.org/wiki/疎行列#圧縮行格納方式_(CRS)・Compressed_Sparse_Row_(CSR))などを参考にしました。

上記ソースコード内の monolish と書いてるところが今回のmonolish関係のコードです。サンプルコード通りにprecondを行うと答えが合わなかったのでコメントアウトしてあります。
`send`/`recv`はGPU版で必要(GPUにデータを送る/GPUからデータを返す)で、CPU版では必要はないです。書いてても動きます。

```bash
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は転送時間を含みます。

![Screen Shot 2021-04-10 at 19.40.52.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/328786/5735ff0f-acc4-3ee8-7c54-a1addbcef057.png)

ちゃんとGPUが使えてて、N=22ではGPUを使って3倍高速化できてそうです。
N=16などの小さいところでは、GPUとの通信のせいでCPUよりも計算が遅くなっていそうです。

N=24はGPUのメモリが足りませんでした。。。
メモリたくさんのGPU欲しい。。。


monolishのインストールに半年で4回計24時間ぐらい失敗し続けたので、この記事が多少の助けになると幸いです。