LoginSignup
11
6

More than 1 year has passed since last update.

物理学アドベントカレンダーAdvent Calendar 2021Day 8
@fockl

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

Last updated at Posted at 2021-12-07

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

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

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

monolishってなぁに

私にとっては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++にしたものがこちらになります。

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;
}

量子スピン模型のハミルトニアンを数値対角化しようの時とは違い、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は転送時間を含みます。

Screen Shot 2021-04-10 at 19.40.52.png

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

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

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

11
6
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
11
6