페이지

2014년 12월 29일 월요일

OpenMP를 이용한 멀티 쓰레드 프로그래밍 HOWTO #2

C언어: OpenMP를 이용한 멀티 쓰레드 프로그래밍 HOWTO #2 




목차
0. 시작하기에 앞서
1. OpenMP란?
2. OpenMP의 시작
3. OpenMP를 이용한 work-sharing 모델
4. Loop construct
  4.a Loop consturct 기본 형태
  4.b 병렬 처리시 주의점
  4.c 스케줄링

* 연습문제 (멀티 쓰레딩과 reentrant 함수의 관계에 대한 문제)
5. Section construct
* 연습문제 (섹션별로 task를 할당하는 방법)
6. Single Construct
7. Master Construct
8. Barrier
  8.a Implicit barrier
  8.b Explicit barrier
9. critical & atomic construct
  9.a critical construct
  9.b atomic construct
10. OpenMP API
11. OpenMP에서 제공하는 환경변수 및 전처리기
12. OpenMP vs pthread (POSIX thread)
13. 마치면서
 


* 연습문제 (멀티 쓰레딩과 reentrant 함수의 관계에 대한 문제)



이번에는 pi를 몬테 카를로 시뮬레이션을 이용해서 구해보도록 하겠다.
각변의 길이가 1인 정사각형이 있다. 그리고 정사각형 안에 반지름(r)이 1인 원호를 그리도록 하자.
그러면 중점으로부터 원호까지의 최단 길이는 무조건 1이 된다.

그러면 이제 정사각형 안에 임의의 점(x,y)좌표를 찍어서 중점에서 선분을 연결하고 아랫변까지 선분을 내리면 직각삼각형이 된다.
그리고 이 직각삼각형의 길이는 x, 높이는 y가 된다.
피타고라스의 정리에 의해 (x의 제곱)+(y의 제곱)=(빗변의 제곱)이 되는데, 빗변이 길이가 1보다 작으면
원호 안에 찍힌 점이 되고, 1보다 크면 원호 밖에 정사각형 내부에 찍힌 점이 된다.

원호의 넓이는 반지름 r=1일 때 pi/4이 된다. 따라서 위의 수많은 랜덤 좌표를 찍은 뒤에 (원호 내부의 점의 개수)/(전체 랜덤 수)
는 pi/4와 같아질 것이다. 그러면 이제 기본 코딩을 해보자.

#include <stdlib.h>
#include <stdio.h>
#include <time.h>
#include <math.h>
#define LOOP_ITERATION 200000000
int hits;

int main()
{
    int i;
    double x, y, rns;
    time_t t_now;

    printf("Loop iteration = %ld\n", (long)LOOP_ITERATION);
    rns = 1.0 / (double)RAND_MAX;

    t_now = time(0);
    srand((unsigned int)t_now);
    for (i = 0; i<LOOP_ITERATION; i++) {
       x = (double)rand() * rns;
       y = (double)rand() * rns;
       if (x*x + y*y < 1) {
           hits++;
       }
    }
    printf("pi = %f\n", 4 * (double)hits / LOOP_ITERATION);
    return 0;
}

이제 컴파일 후 실행을 해 보겠다. 소스코드의 파일명은 pi_monte.c라고 저장했다. 
$ gcc -o pi_monte pi_monte.c $ time ./pi_monte Loop iteration = 200000000 pi = 3.141588 real 0m8.726s user 0m8.702s sys 0m0.023s
CPU를 한 개만 사용했기 때문에 real과 user+sys 시간이 같게 나온다. 이제 여기에 OpenMP를 적용해보도록 하자. 중복되는 코드는 모두 생략하고 중간에 loop부분만 적도록 하겠다. 
...생략...
#pragma omp parallel for private(x,y) reduction(+:hits)
   for (i = 0; i<LOOP_ITERATION; i++) {
       x = (double)rand() * rns;
       y = (double)rand() * rns;
       if (x*x + y*y < 1) {
          hits++;
       }
   }
...생략...

자 이제 수정된 소스코드를 컴파일하고 다시 실행해보겠다. 
 
$ gcc -o pi_monte_omp -fopenmp pi_monte.c $ time ./pi_monte_omp Loop iteration = 200000000 pi = 3.141513 real 0m54.909s user 0m54.173s sys 0m53.115s

OpenMP버전의 실행 시간이 엄청나게 증가한 것을 볼 수 있다.
이는 뭔가 문제가 발생한 것이다. profiler가 있으면 cache miss가 많아진 것을 추적할 수 있다. 그러면 왜 cache miss가 많아졌을까?

이는 멀티 쓰레드 안전(thread-safe, MT-safe)한 함수를 사용하지 않았기 때문이다.
단도직입으로 원인은 rand()함수이다. rand()함수는 내부에 static 변수(BSS영역)를 사용하기 때문에 lock 없이 사용하면
오염된 공간을 쓸 수 있다. 그렇다고 lock으로 보호하는 것은 성능상으로 좋지 못하다. 따라서 MT-safe한 랜덤생성 함수로
대체해야 한다. 마침 rand()의 reentrant 버전인 rand_r()이 있으므로 이를 사용하도록 바꾸면 된다.

* reentrant에 대한 이해가 필요하다면...http://sunyzero.egloos.com/4188321 글을 읽도록 한다.

그러면 이제 수정된 버전을 보도록 하겠다. 

#define _REENTRANT
#include <stdlib.h>
#include <stdio.h>
#include <time.h>
#include <math.h>
#ifdef LOOP_ITERATION
#define LOOP_ITERATION 200000000
#endif
int hits;
int main()
{
    unsigned int state;
    int i;
    double x, y, rns;
    printf("Loop iteration = %ld\n", (long)LOOP_ITERATION);
    rns = 1.0 / (double)RAND_MAX;
    state = (unsigned int)time(0);
#pragma omp parallel for firstprivate(state) private(x, y) reduction(+:hits)
    for (i = 0; i<LOOP_ITERATION; i++) {
       x = (double)rand_r(&state) * rns;
       y = (double)rand_r(&state) * rns;
#if LOOP_ITERATION < 100
       printf("THR[%d:%d] x,y/state = %f,%f/%u\n", omp_get_thread_num(), i, x, y, state);
#endif
       if (x*x + y*y < 1) {
           hits++;
       }
    }
    printf("pi = %f (hits = %d)\n", (double)hits / LOOP_ITERATION * 4, hits);
    return 0;
}

중요한 변화를 확인하기 위해 LOOP_ITERATION이 100 이하인 경우에는 루프를 돌때마다 rand_r()로 얻어진 x, y, state의 값을 출력하도록 디버깅 코드를 심어놨다. 그러면 LOOP_ITERATION을 20으로 줄여서 컴파일&실행으로 디버깅 해보자.

$ gcc -o pi_monte_omp_logging -DLOOP_ITERATION=20 -fopenmp pi_monte.c $ ./pi_monte_omp_logging Loop iteration = 20 THR[0:0] x,y/state = 0.739258,0.373189/1534713604 THR[0:1] x,y/state = 0.741886,0.152093/434332526 THR[0:2] x,y/state = 0.399359,0.675170/1402811528 THR[0:3] x,y/state = 0.140991,0.510551/3217124562 THR[0:4] x,y/state = 0.117010,0.018486/2612351692 THR[0:5] x,y/state = 0.151627,0.286311/535820534 THR[0:6] x,y/state = 0.412198,0.450266/361487824 THR[0:7] x,y/state = 0.126866,0.324974/1313573850 THR[0:8] x,y/state = 0.923697,0.862684/1937324436 THR[0:9] x,y/state = 0.828232,0.810037/3167316350 THR[1:10] x,y/state = 0.739258,0.373189/1534713604 THR[1:11] x,y/state = 0.741886,0.152093/434332526 THR[1:12] x,y/state = 0.399359,0.675170/1402811528 THR[1:13] x,y/state = 0.140991,0.510551/3217124562 THR[1:14] x,y/state = 0.117010,0.018486/2612351692 THR[1:15] x,y/state = 0.151627,0.286311/535820534 THR[1:16] x,y/state = 0.412198,0.450266/361487824 THR[1:17] x,y/state = 0.126866,0.324974/1313573850 THR[1:18] x,y/state = 0.923697,0.862684/1937324436 THR[1:19] x,y/state = 0.828232,0.810037/3167316350 pi = 3.200000 (hits = 16)
iteration이 20번이므로 pi의 정확도는 일단 포기하자. 지금은 pi의 결과를 보려는 것이 아니라, 각 쓰레드의 x, y값을 확인하기 위함이다. 0번째 쓰레드의 첫번째 데이터와 1번째 쓰레드의 첫번째 데이터를 비교하니 뭔가 이상하다. 이해를 돕기 위해 두개만 따로 떼어서 보도록 하자.

THR[0:0] x,y/state = 0.739258,0.373189/1534713604
THR[1:10] x,y/state = 0.739258,0.373189/1534713604

둘이 놀랍도록 일치하지 않은가? 왜 이런 문제가 발생할까? 이는 rand_r()함수에 사용하는 seed값 변수인 state의 초기값이 모든 쓰레드들에 대해서 같기 때문에 발생하는 문제다. 그렇다면 각 쓰레드가 seed를 다르게 가져가도록 소스코드를 수정해야 한다.

#define _REENTRANT
#include <stdlib.h>
#include <stdio.h>
#include <time.h>
#include <math.h>
#ifdef LOOP_ITERATION
#define LOOP_ITERATION 200000000
#endif
int hits;
int main()
{
    unsigned int state1, state2;
    int i;
    double x, y, rns;
    printf("Loop iteration = %ld\n", (long)LOOP_ITERATION);
    rns = 1.0 / (double)RAND_MAX;
    state1 = (unsigned int)times(NULL);
#pragma omp parallel private(x, y, state2) reduction(+:hits) shared(state1)
    {
#pragma omp critical
       state2 = rand_r(&state1);
       printf("THR[%d] state1/state2 = %u/%u\n", omp_get_thread_num(), state1, state2);
#pragma omp for
       for (i = 0; i<LOOP_ITERATION; i++) {
           x = (double)rand_r(&state2) * rns;
           y = (double)rand_r(&state2) * rns;
#if LOOP_ITERATION < 100
           printf("THR[%d:%d] x,y/state,hits = %f,%f/%u,%d\n", omp_get_thread_num(), i, x, y, state2, hits);
#endif
           if (x*x + y*y < 1) {
              hits++;
           }
       }
    }
    printf("hits(%d), pi = %f\n", hits, (double)hits / LOOP_ITERATION * 4);
    return 0;
}

#pragma omp critical 지시어는 아래의 블럭공간을 락으로 보호해준다. 즉 state1으로부터 state2라는 seed값을 생성하는 코드 부분은 각 쓰레드들이 직렬적으로 실행하게 되므로 각각 다른 seed (=state2)를 가질 수 있게 해준다.

그러면 이제 수정된 버전의 확인을 위해 LOOP_ITERATION을 20번으로 지정하고 컴파일하여 다시 실행해보자. 
$ gcc -o pi_monte_omp_logging -DLOOP_ITERATION=20 -fopenmp pi_monte.c $ ./pi_monte_omp_logging Loop iteration = 20 THR[0] state1/state2 = 2209777173/920556694 THR[0:0] x,y/state,hits = 0.007502,0.800638/402955376,0 THR[0:1] x,y/state,hits = 0.031584,0.366846/1879397754,1 THR[0:2] x,y/state,hits = 0.084310,0.126795/2318644788,2 THR[0:3] x,y/state,hits = 0.219188,0.054924/839319838,3 THR[0:4] x,y/state,hits = 0.315477,0.640078/2081973432,4 THR[0:5] x,y/state,hits = 0.423693,0.287051/4130577282,5 THR[0:6] x,y/state,hits = 0.738140,0.376765/2380434428,6 THR[0:7] x,y/state,hits = 0.127375,0.723816/3248220326,7 THR[0:8] x,y/state,hits = 0.832496,0.256793/3481063424,8 THR[0:9] x,y/state,hits = 0.939841,0.108509/3174395018,9 THR[1] state1/state2 = 2209777173/754906038 THR[1:10] x,y/state,hits = 0.446803,0.061819/2639594128,0 THR[1:11] x,y/state,hits = 0.919237,0.585474/2439875226,1 THR[1:12] x,y/state,hits = 0.122634,0.254468/2048737876,1 THR[1:13] x,y/state,hits = 0.204702,0.526382/3058641982,2 THR[1:14] x,y/state,hits = 0.642845,0.756832/477247192,3 THR[1:15] x,y/state,hits = 0.550340,0.531024/1736039586,4 THR[1:16] x,y/state,hits = 0.926036,0.217729/4074357788,5 THR[1:17] x,y/state,hits = 0.266219,0.975137/2510577606,6 THR[1:18] x,y/state,hits = 0.663719,0.557898/2880736800,6 THR[1:19] x,y/state,hits = 0.491536,0.580296/2760522154,7 hits(18), pi = 3.600000
실행 결과를 보면 이제 x, y값이 각 쓰레드마다 달라지는 것을 볼 수 있다.

그러면 이제 기능상의 문제는 해결되었으니 원래대로 200,000,000번의 시행 횟수로 실행시간을 비교해보자. 
 
$ gcc -o pi_monte_omp -fopenmp pi_monte.c $ time ./pi_monte_omp Loop iteration = 200000000 THR[0] state1/state2 = 319398670/265139008 THR[1] state1/state2 = 319398670/1387545353 hits(157082826), pi = 3.141657 real 0m2.414s user 0m4.737s sys 0m0.013s
 
싱글 쓰레드 버전이 8.6초 걸린 것에 비해 OpenMP로 돌린 버전은 2.4초로 확 줄어들었다. 원래대로라면 듀얼코어니 4.3초정도가 나와야 하겠지만, rand_r()자체가 rand()보다 가볍기 때문에 그로 인해서 속도가 더 빨라졌다. 즉 싱글 쓰레드 버전이라고 해도 rand()보다는 rand_r()을 사용하면 더 빠르다라는 교훈도 덤으로 알려드리는 문제였다.

5. Section constructsection construct는 task level parallelism에 사용하며, 각각의 작업들이 서로 관련이 없는 경우에 사용한다.
(task 레벨의 병렬화이므로 divide and conquer 형태의 문제해결에도 적용가능하다.)

그림에서 볼 수 있듯이 오렌지색의 섹션과 블루색의 섹션이 각각 독립적으로 작동하도록 구성할 수 있다.
단 각각의 섹션은 누가 먼저 종료하든지 #pragram omp sections 블록의 끝에는 implicit barrier가 있으므로 대기하게 된다.

참고로 병렬구간내에 섹션구간만 존재하는 경우라면 #pragma omp parallel sections로 구문을 합칠 수 있다.

그러면 이번에는 앞에서 loop construct때 연습했던 2가지 pi 구하는 방법(numerical integration, monte carlo simulation)을
sections를 이용해서 동시에 작동시키도록 바꿔보자.


* 연습문제 (섹션별로 task를 할당하는 방법)




앞에서 Single threaded 버전으로 만들었던 소스코드들을 연달아 붙여두면 된다.
우선 Numerical Integration 방법의 소스코드를 다시 보자. 

#include <stdio.h>
#include <stdlib.h>
int num_steps = 1000000000; /* 10억번 : 너무 많으면 조금 줄이시길... */
int main()
{
    int i;
    double x, step, sum = 0.0;
    step = 1.0 / (double)num_steps;
    for (i = 0; i<num_steps; i++) {
       x = (i + 0.5) * step;
       sum += 4.0 / (1.0 + x*x);
    }
    printf("PI = %.8f (sum = %.8f)\n", step*sum, sum);
    return EXIT_SUCCESS;
}

이번에는 Monte Carlo simulation방법을 다시 보겠다. 
#include <stdlib.h>
#include <stdio.h>
#include <time.h>
#include <math.h>
#define LOOP_ITERATION 200000000
int hits;
int main()
{
    int i;
    double x, y, rns;
    time_t t_now;
    printf("Loop iteration = %ld\n", (long)LOOP_ITERATION);
    rns = 1.0 / (double)RAND_MAX;
    t_now = time(0);
    srand((unsigned int)t_now);
    for (i = 0; i<LOOP_ITERATION; i++) {
       x = (double)rand() * rns;
       y = (double)rand() * rns;
       if (x*x + y*y < 1) {
           hits++;
       }
    }
    printf("pi = %f\n", 4 * (double)hits / LOOP_ITERATION);
    return 0;
}

이제 이 2개의 소스코드를 합쳐야 한다. 주의할 점은 두번째 Monte Carlo simulation에서는 rand()함수 대신에 rand_r()을 사용하는 것을 잊지 말아야 한다. 항상 Non-multi-threaded model을 Multi-threaded model로 바꿀때는 MT-safe function이나 reentrant function으로 골라 써야 한다는 점이다.

그리고 false-sharing문제를 피하기 위해서 동일 캐시 라인에 올라갈 수 있는 전역변수나 힙 메모리를 사용하면 안된다는 점이다. 왠만하면 쓰레드 로컬 변수에서 대부분 해결하도록 하자.

그러면 둘을 #pragma omp sections로 합친 코드를 보겠다. 보기 좋게 출력부분의 수정을 했으나 기본 코드는 같다.
#define _XOPEN_SOURCE   600
#include <stdlib.h>
#include <unistd.h>
#include <stdio.h>
#include <time.h>
#include <sys/times.h>
#include <omp.h>
#ifndef LOOP_ITERATION
#define LOOP_ITERATION 200000000
#endif
double  pi[2];
int main()
{
    clock_t     start_m, end_m;
    clock_t     sc_clk_tck = sysconf(_SC_CLK_TCK);
    start_m = times(NULL);
    printf("LOOP ITERATION = %ld\n", (long)LOOP_ITERATION);
#ifdef _OPENMP
    omp_set_num_threads(2);
#endif
#pragma omp parallel
    {
#pragma omp sections
    {
#pragma omp section
    {
       int     i;
       double  x, step, hits;
       float   t_elapsed;
       clock_t start, end;
       hits = 0.0;
#ifdef _OPENMP
       printf("[SECTION] integration method by thread(%d)\n", omp_get_thread_num());
#else
       printf("[SECTION] start integration method\n");
#endif
       start = times(NULL); /* get clock tick */
       step = 1.0 / (double)LOOP_ITERATION;
       for (i = 0; i<LOOP_ITERATION; i++) {
           x = (i + 0.5) * step;
           hits += 1.0 / (1.0 + x*x);
       }
       pi[0] = 4 * hits * step;
       end = times(NULL);
       t_elapsed = (float)(end - start) / sc_clk_tck;
#ifdef _OPENMP
       printf("[SECTION] end integration method by thread(%d):elapsed time(%.02f sec)\n",
           omp_get_thread_num(), t_elapsed);
#else
       printf("[SECTION] end integration method: elapsed time(%.02f sec)\n"t_elapsed);
#endif
    }
#pragma omp section
    {
       int     i, state, hits;
       double  x, y, rns;
       float   t_elapsed;
       clock_t start, end;
#ifdef _OPENMP
       printf("[SECTION] monte carlo method by thread(%d)\n", omp_get_thread_num());
#else
       printf("[SECTION] start monte carlo method\n");
#endif
       start = times(NULL); /* get clock tick */
       state = time(0);
       rns = 1.0 / (double)RAND_MAX;
       hits = 0;
       for (i = 0; i<LOOP_ITERATION; i++) {
           x = (double)rand_r((unsigned int *)&state) * rns;
           y = (double)rand_r((unsigned int *)&state) * rns;
           if (x*x + y*y < 1) {
              hits++;
           }
       }
       pi[1] = (hits / LOOP_ITERATION) * 4;
       end = times(NULL);
       t_elapsed = (float)(end - start) / sc_clk_tck;
#ifdef _OPENMP
       printf("[SECTION] end monte carlo method by thread(%d): elapsed time(%.02f sec)\n",
           omp_get_thread_num(), t_elapsed);
#else
       printf("[SECTION] end monte carlo method: elapsed time(%.02f sec)\n", t_elapsed);
#endif
    }
    }
    }
    end_m = times(NULL);
    printf("integration PI = %.8f\n", pi[0]);
    printf("monte carlo PI = %.8f\n", pi[1]);
    printf("* Total elapsed time(%.02f sec)\n", (double)(end_m - start_m) / sc_clk_tck);
    return EXIT_SUCCESS;
}

이제 실행을 해보자. 

$ gcc -o pi_sections_omp -fopenmp pi_sections.c $ time ./pi_sections_omp LOOP ITERATION = 200000000                                    [SECTION] integration method by thread(0)                     [SECTION] monte carlo method by thread(1)                     [SECTION] end integration method by thread(0):elapsed time(1.79 sec) [SECTION] end monte carlo method by thread(1): elapsed time(5.05 sec) integration PI = 3.14159265                                          monte carlo PI = 3.14161276                                          * Total elapsed time(5.05 sec)                                       real    0m5.052s user    0m6.858s sys     0m0.016s
당연히 Numerical Integration이 더 빠르기 때문에 먼저 끝난다. 하지만 implicit barrier가 있기 때문에 대기하게 된다. Monte Carlo simulation은 더 오래 걸리기 때문에 전체 수행 시간은 Monte Carlo simulation이 끝나는 시간에 종료한다.

* 오늘은 여기까지... 나머지는 다음 장에. (이것도 정말 힘들군요... 본문의 반말은 이해해주세요.)


댓글 없음:

댓글 쓰기