続：2次元イジングモデル

前回の2次元イジングモデルのコードをCで書きなおしてみました。unixの標準コマンドtimeを使って、openmpありとなしの実行時間を計算してみます。2次元の行列サイズは1024x1024、ループ回数は1000です。コード中のMT.hは乱数を生成するためにメルセンヌツイスターを使用しているからです。

環境

• OS:Ubuntu 14.04
• CPU:Core i7 5820K
• Mem:CT2K8G4DFD8213
• Mother:ASRock X99 Extreme4
• コンパイラオプション:icc -openmp -O3

結果

openmp	time
Yes(10cores)	4.017s
No	10.963s

コード

#include <stdio.h>
#include <stdlib.h>
#include "MT.h"
#include <omp.h>
#include <math.h>

float calcE(int *m,int row,int col){
    int i,j,nx,ny;
    float J = 1;
    float E = 0;
    #pragma omp parallel num_threads(10)
    {
    #pragma omp for private(i,j),reduction(+:E)
    for(i=0;i<col;i++){
        for(j=0;j<row;j++){
            nx = j;
            ny = (i+1)%col;
            E += -J * m[i*col+j]*m[ny*col+nx];

            nx = (j+1)%row;
            ny = i;
            E += -J * m[i*col+j]*m[ny*col+nx];
            
            nx = j;
            ny = ((i-1)%col != -1) ? i : col;
            E += -J * m[i*col+j]*m[ny*col+nx];

            nx = ((j-1)%row != -1) ? j : row;
            ny = i;
            E += -J * m[i*col+j]*m[ny*col+nx];
        }
    }
    #pragma omp barrier
    }
    return E;
}

int copy(int *dst,int *src,int size){
    int i;
    for(i=0;i<size;i++)
        dst[i]=src[i];
    return(0);
}

int init(int *m,int size){
    int i;
    int sign[] = {-1,1};

    for(i=0;i<size;i++)
        m[i] = sign[genrand_int32()%2];
    return(0);
}

void MC(int *m,int row,int col,int nstep){
    int i;
    int mp;
    int size = row*col;
    int *tmpm = (int *) malloc(sizeof(int)*size);
    float Epre,Eafter;
    float prob;
    float T = 0.5;
    float invT = 1/T;
    
    for(i=0;i<nstep;i++){
        mp = genrand_int32()%size;
        copy(tmpm,m,size);
        tmpm[mp] = -tmpm[mp];
        Epre = calcE(m,row,col);
        Eafter = calcE(tmpm,row,col);

        if (Eafter < Epre)
            copy(m,tmpm,size);
        else{
            prob = genrand_real1();
            if(prob < exp(-(Eafter-Epre)*invT))
                copy(m,tmpm,size);
        }
        printf("%lf\n",Eafter);
    }
}

int main(void){
    int bsize = 10;
    int row = 1<<bsize;
    int col = 1<<bsize;
    int size = row * col;
    int nstep = 1000;
    int seed = 1;
    int *mat;    

    init_genrand(seed);

    mat = (int *) malloc(sizeof(int)*size);
    init(mat,size);

    MC(mat,row,col,nstep);
    free(mat);
    return(0);
}

2015/11/12:追記

float calcE(int *m,int row,int col){
    int i,j,nx,ny;
    float J = 1;
    float E = 0;
    #pragma omp parallel num_threads(10)
    {
    #pragma omp for private(i,j),reduction(+:E)
    for(i=0;i<col;i++){
        for(j=0;j<row;j++){
            nx = j;
            ny = (i+1)%col;
            E += -J * m[i*col+j]*m[ny*col+nx];

            nx = (j+1)%row;
            ny = i;
            E += -J * m[i*col+j]*m[ny*col+nx];
            
            nx = j;
            ny = ((i-1)%col != -1) ? i : col;
            E += -J * m[i*col+j]*m[ny*col+nx];

            nx = ((j-1)%row != -1) ? j : row;
            ny = i;
            E += -J * m[i*col+j]*m[ny*col+nx];
        }
    }
    #pragma omp barrier
    }
    return E;
}

は正しくは

float calcE(int *m,int row,int col){
    int i,j,nx,ny;
    float J = 1;
    float E = 0;
    #pragma omp parallel num_threads(10)
    #pragma omp for private(i,j),reduction(+:E)
    for(i=0;i<row;i++){
        for(j=0;j<col;j++){
            nx = i;
            ny = (j+1)%col;
            E += -J * m[i*row+j]*m[nx*row+ny];

            nx = (i+1)%row;
            ny = j;
            E += -J * m[i*row+j]*m[nx*row+ny];

            nx = i;
            ny = ((j-1)%col != -1) ? j-1 : col-1;
            E += -J * m[i*row+j]*m[nx*row+ny];
            
            nx = ((i-1)%row != -1) ? i-1 : row-1;
            ny = j;
            E += -J * m[i*row+j]*m[nx*row+ny];
        }
    }
    #pragma omp barrier

    return E;
}