続*6：2次元イジングモデル
現在、CUDAで２次元イジングモデルのプログラムを書いています。下記のgifは、プログラミングミスの結果ですが、面白かったので載せておきます。なんだかライフゲームみたいです。
12/21追記
CUDAで2次元イジングモデルのプログラムを書きました。いろいろと直すべき点はありますが、まずは実行できるモデルを作成しました。外部地場Hは0のモデルです。まずは、2024*2024の系のgifから
結果
| コンパイラオプション(GTX 760) | nvcc -O3 |
|---|---|
| 環境(3072*3072サイズ):Time | 23.919s |
コードの説明
ヘッダファイル
#include <iostream> #include <cuda.h> #include <cuda_runtime.h> #include <curand.h> #include <curand_kernel.h> #include <stdio.h>
今回は、C++でプログラミングしました。ここでのポイントはcurand.hとcurand_kernel.hです。今回は、curand_kernel.hしか使用していませんが、教育的な意味でcurand.hもインクルードしています(curand.hの存在があることを思い出すため)。curand.hとcurand_kernel.hは、どちらも乱数列を生成するためのモジュール郡です。その違いは、curand.hはホストから呼び出し可能な関数群が格納されており、curand_kernel.hはデバイスからのみ呼び出し可能な関数群を格納しています。つまり、乱数列のみが欲しい場合は、cudaプログラミングをする必要がなく、curand.hを使って高速に乱数列を生み出すことができるのです。
その他の宣言
using namespace std; #define D_CHECK(err){\ if(err !=cudaSuccess)\ printf("error %s,at %d\n",cudaGetErrorString(err),__LINE__); \ } #define ROW 2048 #define COL 2048 __global__ void g_rand_init(int size,curandState *states); __global__ void g_calc_energy(int J,int *S,int *E,int row,int col,curandState *states); __global__ void g_S_init(int *S,int size,curandState *states); __global__ void g_simulate(int J,float invT,int *S,int *E,int row,int col,curandState *states,int flag);
D_CHECKは、cudaの関数を呼び出した時のエラーチェック用の簡単な宣言です。まず、ROWとCOLでイジングの数を指定しています。その次に、CUDAのカーネル関数を宣言しています。注意するべき点は、カーネル関数はメンバ関数にすることができません。そのため、外部で宣言しています。
クラスの宣言
template <typename T> class ising{ public: // host values int row; int col; int size; T J; // 交換相互作用エネルギー int block_x,block_y; int thread_x,thread_y; dim3 grid,block; int nthreads; int sumE; float invT; FILE *fp; T *S; // スピンの情報 T *E; // 各スピンのエネルギー T *hE; // エネルギー計算用 // device values cudaError_t err; curandState *states; //乱数列生成用 T *dS; // デバイス側のスピン T *dE; // デバイス側の各スピンのエネルギー //functions ising(int row_,int col_,T J_,float temperature); ~ising(); void devInit(int bx,int by,int tx,int ty); // デバイス側の初期化 void toFile(); // ファイルの書き出し void checkEnergy(); // 系のエネルギーのチェック void run(); //シミュレーションのスタート void devEnd(); // デバイスの変数の解放 };
templateを使用しているのは、今後、XYモデルやハイゼンベルグモデルなんかの実装も考えているからです。ここでcurandStateで宣言しているstatesが各スレッドでの乱数列の情報を持ちます。
コンストラクタとデコンストラクタ
template <typename T> ising<T>::ising(int row_,int col_,T J_,float temperature_){ row = row_; col = col_; size = row_ * col_; J = J_; invT = 1.0/temperature_; sumE = 0; S = (T *) malloc(sizeof(T) * size); E = (T *) malloc(sizeof(T) * size); hE = (T *) malloc(sizeof(T) * size); if((fp = fopen("result.dat","w")) == NULL){ printf("Error at opening file.\n"); exit(0); } fprintf(fp,"#%9d%9d\n",row,col); for(int i=0;i<size;i++) S[i] = 0; } template <typename T> ising<T>::~ising(){ free(S); free(E); free(hE); fclose(fp); }
デバイスの初期化と解放
template <typename T> void ising<T>::devInit(int bx,int by,int tx,int ty){ nthreads = bx * by * tx * ty; dim3 grid_(bx,by),block_(tx,ty); grid = grid_; block = block_; printf("grid=[%d,%d],block[%d,%d]\n",grid.x,grid.y,block.x,block.y); cudaMalloc((curandState **)&states,sizeof(curandState)*nthreads); cudaMalloc((T **)&dS,sizeof(T)*size); cudaMalloc((T **)&dE,sizeof(T)*size); g_rand_init<<<grid,block>>>(size,states); g_S_init<<<grid,block>>>(dS,size,states); D_CHECK(cudaMemcpy(S,dS,sizeof(int)*size,cudaMemcpyDeviceToHost)); } template <typename T> void ising<T>::devEnd(){ cudaFree(states); cudaFree(dS); cudaFree(dE); cudaDeviceReset(); }
スピンの初期化と乱数の初期化
__global__ void g_rand_init(int size,curandState *states){ int x = blockIdx.x * blockDim.x + threadIdx.x; int y = blockIdx.y * blockDim.y + threadIdx.y; int idx = x + y*gridDim.x*blockDim.x; if(idx < size){ curand_init(idx,0,0,&states[idx]); } } __global__ void g_S_init(int *S,int size,curandState *states){ int x = blockIdx.x * blockDim.x + threadIdx.x; int y = blockIdx.y * blockDim.y + threadIdx.y; int idx = x + y*gridDim.x*blockDim.x; if(idx < size){ if(curand_uniform(&states[idx]) <= 0.5) S[idx] = -1; else S[idx] = 1; } }
シミュレーションの中心部分
template <typename T> void ising<T>::run(){ int flag; flag = 0; g_simulate<<<grid,block>>>(J,invT,dS,dE,row,col,states,flag); flag = 1; g_simulate<<<grid,block>>>(J,invT,dS,dE,row,col,states,flag); } __global__ void g_simulate(int J,float invT,int *S,int *E,int row,int col,curandState *states,int flag){ int x = blockIdx.x * blockDim.x + threadIdx.x; int y = blockIdx.y * blockDim.y + threadIdx.y; int idx = x + y*gridDim.x*blockDim.x; int nx; int ny; float prob; if( x < row && y < col){ if((x+y)% 2 == flag){ S[x*col+y] = -S[x*col+y]; E[x*col+y] = 0; nx = x-1; ny = y; if(nx < 0) nx = col-1; E[x*col+y] += -J*S[x*col+y]*S[nx*col+ny]; nx = x; ny = y-1; if(ny < 0) ny = row-1; E[x*col+y] += -J*S[x*col+y]*S[nx*col+ny]; nx = x+1; ny = y; if(nx >= row) nx = 0; E[x*col+y] += -J*S[x*col+y]*S[nx*col+ny]; nx = x; ny = y+1; if(ny >= col) ny = 0; E[x*col+y] += -J*S[x*col+y]*S[nx*col+ny]; if( E[x*col+y]>=0){ prob = curand_uniform(&states[idx]); if(prob > exp(-2.0*E[x*col+y]*invT)){ S[x*col+y] = -S[x*col+y]; } } } } __syncthreads(); }
メイン関数
int main(void){ int row = ROW; int col = COL; int J = 1; ising<int> ising2d(row,col,J,0.5); ising2d.devInit(64,64,32,32); ising2d.toFile(); for(int i=0;i<100;i++){ ising2d.run(); //ising2d.checkEnergy(); //ising2d.toFile(); } ising2d.devEnd(); return 0; }
