[OpenMP] 사다리꼴 공식 (trapezoidal rule)

Refenences

• An Introduction to Parallel Programming

Contents

• 사다리꼴 공식 (trapezoidal rule)
• ciritical 디렉티브
• 변수의 범위(scope)
• reduction 디렉티브
• parallel for 디렉티브

---

2021.11.09 - [프로그래밍/병렬프로그래밍] - [MPI] 사다리꼴 공식 (trapezoidal rule) 병렬화 - 1

이번 포스팅에서는 MPI에서 살펴봤던 사다리꼴 공식을 OpenMP로 어떻게 구현할 수 있는지 다양한 버전으로 알아보겠습니다. 

 

사다리꼴 공식은 x축과 x=a와 b의 f(x) 그래프 사이의 면적을 구하는 공식으로, x축을 n개의 서브 인터벌로 나누고 사다리꼴 넓이 공식을 통해 각 서브 인터벌의 면적을 추정하는 것입니다.

이전 포스팅에서 자세하게 설명했었는데, 각각의 서브 인터벌은 같은 길이 \(h = \frac{b-a}/n\)이며, \(x_i = a + ih, i = 0, 1, \cdots, n\)이라고 했을 때, 그래프 아래의 넓이는 

\[h[\frac{f(x_0)}{2} + f(x_1) + f(x_2) + \cdots + f(x_{n-1}) + \frac{f(x_n)}{2}]\]

으로 구할 수 있습니다.

 

위 수식을 시리얼 코드로 작성하면 다음과 같습니다.

시리얼로 작성된 코드는 다음을 참조하시면 됩니다. 따로 살펴보지는 않겠습니다.

https://github.com/junstar92/parallel_programming_study/blob/master/mpi/01_serial_trap.c

GitHub - junstar92/parallel_programming_study

---

첫 번째 버전

우선 이 사다리꼴 공식을 병렬 프로그램으로 디자인하면, 두 종류의 태스크를 정의할 수 있습니다.

하나는 각 인터벌의 사다리꼴 넓이를 계산하는 것이고, 다른 하나는 계산한 넓이를 글로벌 변수에 더하는 것입니다. 또한, 첫 번째 태스크에서는 태스크 간 통신을 할 필요는 없지만, 두 번째 태스크에서는 첫 번째 태스크의 결과를 얻기 위해 통신이 필요합니다. 그리고, 코어 수보다 사다리꼴의 수가 더 많다고 가정하고, 각 스레드에 사다리꼴 블록을 계속 할당하여 전체 사다리꼴 넓이를 계산합니다.

그러나, 두 번째 태스크처럼 각 스레드들의 결과를 누적시켜야하기 때문에 모든 작업이 끝난 것은 아닙니다. 확실한 솔루션은 모든 스레드들의 결과를 하나의 공유 변수에 누적시키는 것입니다. 그렇다면, 각 스레드는 다음과 같은 문장을 실행하면 됩니다.

global_result += my_result;

잘 아시겠지만, 위 코드는 오류를 발생시키게 됩니다. 두 개 혹은 그 이상의 스레드가 동시에 이 문장을 실행한다면 결과는 예측할 수 없습니다. 

이는 크리티컬 섹션이 보호되지 않아서 발생하고, 크리티컬 섹션의 자세한 내용은 아래 포스트에서 간단하게 설명하고 있습니다 !

2021.11.17 - [프로그래밍/병렬프로그래밍] - [pthread] Critical Section

[pthread] Critical Section

 

그렇기 때문에 한 번에 하나의 스레드가 

global_result += my_result;

를 실행할 수 있도록 하는 방법이 필요하며, 다른 스레드들은 이 스레드가 위 코드 수행을 완료할 때까지는 해당 코드를 실행하지 않도록 해야합니다. Pthreads에서는 뮤텍스나 세마포어를 사용했습니다.

 

OpenMP에서는 critical 디렉티브를 사용합니다.

#pragma omp critical
global_result += my_result;

이 디렉티브는 시스템으로 하여금 스레드들이 구조화된 코드 블록을 상호 배제(mutually exclusive)로 실행하도록 컴파일러에게 지시합니다. 즉, 오직 한번에 하나의 스레드만 디렉티브 뒤에 따라오는 코드 블록을 실행할 수 있습니다.

 

critical 디렉티브를 사용하여 사다리꼴 공식 프로그램의 코드는 다음과 같습니다. 전체 코드는 아래 링크를 참조해주세요. 

https://github.com/junstar92/parallel_programming_study/blob/master/OpenMP/02_omp_trap1.c

GitHub - junstar92/parallel_programming_study

 

#include <stdio.h>
#include <stdlib.h>
#include <omp.h>

void Usage(char* prog_name);
double f(double x); /* function we're integrating */
void Trap(double a, double b, int n, double* global_result_p);

int main(int argc, char* argv[])
{
    if (argc != 2)
        Usage(argv[0]);

    int thread_count = strtol(argv[1], NULL, 10);
    
    double a, b;
    int n;
    printf("Enter a, b, and n\n");
    scanf("%lf %lf %d", &a, &b, &n);
    if (n % thread_count != 0)
        Usage(argv[0]);

    double global_result = 0.0;
#pragma omp parallel num_threads(thread_count)
    Trap(a, b, n, &global_result);

    printf("With n = %d trapezoids, our estimate\n", n);
    printf("of the integral from %f to %f = %f\n", a, b, global_result);

    return 0;
}

void Trap(double a, double b, int n, double* p_global_result)
{
    double h, local_a, local_b;
    int local_n;
    int my_rank = omp_get_thread_num();
    int thread_count = omp_get_num_threads();

    h = (b-a)/n;
    local_n = n/thread_count;
    local_a = a + my_rank*local_n*h;
    local_b = local_a + local_n*h;

    double my_result = (f(local_a) + f(local_b))/2.0;
    for (int i = 1; i < local_n; i++)
        my_result += f(local_a + i*h);
    my_result = my_result*h;

#pragma omp critical
    *p_global_result += my_result;
}

코드를 살펴보겠습니다.

main 함수 line 23 이전까지는 싱글 스레드이며 스레드의 개수나 입력(a, b, n) 값을 얻어 옵니다. line 24의 parallel 디렉티브는 Trap 함수가 thread_count 개의 스레드를 실행하도록 설정합니다. Trap 함수가 호출된 후에 리턴한 후에는 parallel 디렉티브에 의해서 시작된 스레드들은 모두 종료하고 프로그램은 다시 하나의 스레드로 코드를 계속해서 수행합니다. 

 

Trap 함수에서 각 스레드는 OpenMP 함수를 통해 본인의 rank와 parallel에 의해 시작된 Team 안에 존재하는 스레드의 전체 개수를 얻어 옵니다. 각 스레드는 다음 내용들을 결정합니다.

• 사다리꼴 서브 인터널의 길이 (line 40)
• 각 스레드에 할당된 사다리꼴의 개수 (line 41)
• 각 스레드에서 왼쪽 끝 포인트와 오른쪽 끝 포인트 (line 42-43)
• 계산 결과를 global_result에 반영 (line 45-51)

Trap 함수의 끝에서 각 스레드의 결과를 global_result에 더하면 스레드는 종료합니다.

 

---

변수의 범위

시리얼 프로그래밍에서 변수의 범위는 그 변수가 프로그램의 어느 부분에서 사용되느냐에따라 달라집니다. 예를 들어, C함수의 처음에 선언된 변수는 '함수 영역(function-wide)'(지역 변수)라는 범위를 갖게 되며, 함수 내부에서 액세스할 수 있습니다. 반면에 .c 파일 처음 부분에 선언된 변수는 함수 밖에서도 액세스 가능한 '파일 영역(file-wide)'(전역 변수)이라는 범위를 갖습니다. 파일의 어떤 함수에서도 액세스할 수 있습니다.

OpenMP에서는 parallel 블록에서 스레드들이 변수에 액세스할 수 있는지를 참조합니다. 한 팀에 있는 모든 스레드에 의해서 액세스 가능한 변수는 공유(shared) 범위를 갖게 되며, 하나의 스레드에 의해서 액세스 가능한 변수는 개인(private) 범위를 갖습니다.

 

이전 포스트의 "Hello, OpenMP" 프로그램에서 각 스레드에서 사용된 변수(my_rank, thread_count)는 Hello 함수에서 선언되었고 Hello 함수는 parallel 블록 안에서 호출됩니다. 결론적으로 각 스레드에서 사용되는 변수는 스레드의 (private) 스택에 할당되고, 따라서 이 변수들은 모두 private 범위를 갖습니다. 이는 위의 사다리꼴 공식 예제 코드에서도 적용됩니다. parallel 블록에서 Trap 함수를 호출했기 때문에 각 스레드에 의해 사용되는 Trap 함수 내에 존재하는 변수들은 모두 스레드의 스택에 할당됩니다.

그러나 main 함수에서 선언된 변수들(a, b, c, global_result, thread_count)는 parallel 디렉티브에 의해 시작된 모든 스레드들이 액세스할 수 있습니다. 이 변수들의 기본적인 범위는 shared이기 때문입니다. 이를 묵시적으로 사용하고 있는데, 각 스레드는 Trap 함수 호출할 때, a, b, 그리고 n의 값을 읽어옵니다. 이 호출은 parallel 블록 내에서 일어나고, 이 값들이 각각에 대응하는 매개변수로 복사될 때, a, b, n에 액세스합니다.

Trap 함수에서 p_global_result가 private 변수라고 할 지라도, 이 포인터는 parallel 블록 전에 main 함수에서 선언된 global_result 변수의 주소를 참조하므로, global_result는 parallel 블록 수행 후에 출력하는 결과값을 저장하는데 사용됩니다.

*p_global_result += my_result;

따라서, *p_global_result는 share 범위입니다. 만약 이 변수가 각 스레드별로 private 범위를 갖는 변수라면 ciritical 디렉티브는 필요없습니다. 또한, private 영역이라면 parallel 블록이 완료된 후에 main에서 global_result 값을 결정할 수 없습니다.

 

OpenMP는 변수의 기본 범위를 변경하는 clause를 제공합니다. 다른 디렉티브를 사용하여 변수의 기본 범위를 변경하는 방법을 알아보겠습니다.

---

두 번째 버전

시리얼 프로그램에서 사다리꼴 공식을 적용하여 계산할 때, 그 결과값을 함수의 리턴값으로 전달할 수 있습니다.

double Trap(double a, double b, int n);

그러면 함수는 다음과 같이 호출할 수 있습니다.

global_result = Trap(a, b, n);

 

OpenMP에서도 각 스레드의 로컬 연산을 누적해야하기 때문에, 이와 같이 Local_trap 함수를 정의하여 구현할 수도 있습니다.

double Local_trap(double a, double b, int n)
{
    double h, local_a, local_b;
    int local_n;
    int my_rank = omp_get_thread_num();
    int thread_count = omp_get_num_threads();

    h = (b-a)/n;
    local_n = n/thread_count;
    local_a = a + my_rank*local_n*h;
    local_b = local_a + local_n*h;

    double my_result = (f(local_a) + f(local_b))/2.0;
    for (int i = 1; i < local_n; i++)
        my_result += f(local_a + i*h);
    my_result = my_result*h;

    return my_result;
}

눈에 띄는 특징으로는 이 함수에는 크리티컬 섹션이 없습니다. 각 스레드는 연산의 결과로 my_result 변수의 값을 리턴합니다. 따라서, main 함수는 다음과 같이 수정됩니다.

int main(int argc, char* argv[])
{
    if (argc != 2)
        Usage(argv[0]);

    int thread_count = strtol(argv[1], NULL, 10);
    
    double a, b;
    int n;
    printf("Enter a, b, and n\n");
    scanf("%lf %lf %d", &a, &b, &n);
    if (n % thread_count != 0)
        Usage(argv[0]);
    double start, finish;
    start = omp_get_wtime();
    double global_result = 0.0;
#pragma omp parallel num_threads(thread_count)
    {
#pragma omp critical
        global_result += Local_trap(a, b, n);
    }
    finish = omp_get_wtime();

    printf("With n = %d trapezoids, our estimate\n", n);
    printf("of the integral from %f to %f = %f\n", a, b, global_result);
    printf("Elapsed time = %f seconds\n", finish-start);

    return 0;
}

위 코드에 문제가 하나 있는데, 혹시 눈치 채셨나요?

크리티컬 섹션은 

global_result += Local_trap(double a, double b, int n);

으로 설정되어서, Local_trap의 호출이 한 번에 하나의 스레드에 의해서만 수행됩니다. 이 버전의 런타임을 체크해보면 하나의 스레드를 실행하는 것보다 여러 개의 스레드로 실행하는 것이 더 느릴 수 있습니다.

물론 parallel 블록 안에 private 변수를 선언하고, 다음의 코드처럼 수정하면 크리티컬 섹션을 이동시켜 이러한 문제를 피할 수 있습니다.

int main(int argc, char* argv[])
{
    if (argc != 2)
        Usage(argv[0]);

    int thread_count = strtol(argv[1], NULL, 10);
    
    double a, b;
    int n;
    printf("Enter a, b, and n\n");
    scanf("%lf %lf %d", &a, &b, &n);
    if (n % thread_count != 0)
        Usage(argv[0]);
    double start, finish;
    start = omp_get_wtime();
    double global_result = 0.0;
#pragma omp parallel num_threads(thread_count)
    {
        double my_result = 0.0;
        my_result += Local_trap(a, b, n);
#pragma omp critical
        global_result += my_result;
    }
    finish = omp_get_wtime();

    printf("With n = %d trapezoids, our estimate\n", n);
    printf("of the integral from %f to %f = %f\n", a, b, global_result);
    printf("Elapsed time = %f seconds\n", finish-start);

    return 0;
}

Local_trap의 호출은 크리티컬 섹션 바깥에 있기 때문에 각 스레드들은 동시에 Local_trap을 호출할 수 있습니다.

 

---

세 번째 버전 (Reduction Clause)

OpenMP는 두 번째 버전처럼 Local_trap을 한 번에 하나씩 호출하는 것을 피하기 위한 방법을 제공합니다.

바로 reduction clause인데, global_result를 reduction 변수로 지정하게 되면, 해당 변수는 Thread별로 독립 공간을 유지한 채로 코드를 수행하게 됩니다. 그리고 Thread가 종료되면, Thread별로 계산한 변수 값들을 모아서 한번에 다시 계산합니다.

global_result = 0.0;
#pragma omp prallel num_threads(thread_count) \
	reduction(+: global_result)
global_result += Local_trap(a, b, n);

위와 같이 작성할 수 있습니다. (디렉티브도 백슬래시, '\',를 사용하여 여러 줄로 작성이 가능합니다.)

위 코드는 global_result가 reduction variable이며, reduction( 다음에 나오는 연산자는 reduction operator가 덧셈이라는 것을 나타냅니다. 이 코드가 실행되면, OpenMP는 각 스레드를 위한 private 변수를 생성하며 각 스레드의 결과를 이 private 변수에 저장합니다. 따라서 Local_trap이 병렬적으로 호출될 수 있습니다.

 

reduction clause의 문법은 다음과 같습니다.

reduction(<operator>:<variable list>)

operator는 +, *, -, &, |, ^, &&, || 중의 하나이고, 뺄셈의 경우에는 결합이나 교환 법칙이 성립하지 않아 약간의 문제가 있긴 합니다.

result = 0;
for (int i = 1; i <= 4; i++)
	result -= i;

위 코드를 reduction으로 병렬 수행할 경우 -1-2-3-4 = -10이 되도록 해야합니다. 하지만 OpenMP 표준에서는 이와 같은 결과값을 보장하지는 않습니다.

 

reduction 변수가 float나 double인 경우에도 약간의 문제가 발생할 수 있습니다. 이 결과는 스레드의 숫자에 따라서 약간 달라지는데, 이는 부동소수점 연산에 결합 법칙이 성립하지 않기 때문입니다. 예를 들어, a,b,c가 float형일 때, (a + b) + c는 a + (b + c)와 정확하게 같지 않을 수 있습니다.

 

변수가 reduction clause에 추가될 때 그 변수는 공유 범위를 갖습니다. 하지만 private 변수가 각 스레드에서 생성되고, parallel 블록에서 스레드를 수행할 때에는 private 변수를 사용합니다. 그리고 parallel 블록이 끝나면 private 변수의 값은 공유 변수의 값과 합쳐집니다.

결과적으로 세 번째 버전의 코드는 두 번째 버전과 동일하다고 볼 수 있습니다.

 

reduction의 한 가지 특징은 스레드에 생성되는 private 변수의 초기값은 operator에 따라서 다르게 초기화된다는 점입니다. operator가 +나 -인 경우에는 0으로 초기화되고, 곱셈의 경우에는 1로 초기화됩니다.

 

전체 코드는 아래 링크 참조바랍니다.

https://github.com/junstar92/parallel_programming_study/blob/master/OpenMP/04_omp_trap3.c

GitHub - junstar92/parallel_programming_study

---

네 번째 방법 (parallel for 디렉티브)

사다리꼴 공식의 병렬화 방법으로 OpenMP에서 제공하는 parallel for 디렉티브를 사용할 수 없습니다. 이 디렉티브를 사용하면 다음의 시리얼로 된 프로그램을 간단하게 병렬화할 수 있습니다.

바로 루프전에 디렉티브만 추가해주면 됩니다.

parallel 디렉티브처럼 parallel for 디렉티브도 스레드를 포크(fork)해서 그 이후에 나오는 구조화된 블록을 실행합니다. 그러나 parallel for 디렉티브 다음에 나오는 블록은 반드시 for문이어야 합니다. parallel for 디렉티브를 사용하면 시스템은 루프의 반복을 스레드로 나누어 for 루프를 병렬화합니다.

 

parallel for 디렉티브로 병렬화된 for 루프의 기본적인 분할 방법은 설정한 스레드의 수만큼 스레드를 생성하여 분할합니다. 만약 thread_count개의 스레드를 사용하고 루프가 m번 반복한다면 m/thread_count의 첫 번째 값을 스레드 0에 할당하고, 그 다음 m/thread_count를 스레드 1에 할당하는 방법으로 분할합니다.

 

parallel for을 사용하여 코드를 작성하면 다음과 같습니다.

int main(int argc, char* argv[])
{
    if (argc != 2)
        Usage(argv[0]);

    int thread_count = strtol(argv[1], NULL, 10);
    
    double a, b;
    int n;
    printf("Enter a, b, and n\n");
    scanf("%lf %lf %d", &a, &b, &n);

    double global_result = 0.0;
    global_result = Trap(a, b, n, thread_count);

    printf("With n = %d trapezoids, our estimate\n", n);
    printf("of the integral from %f to %f = %f\n", a, b, global_result);

    return 0;
}

double Trap(double a, double b, int n, int thread_count)
{
    double h, approx;

    h = (b-a)/n;
    approx = (f(a) + f(b))/2.0;
#pragma omp parallel for num_threads(thread_count) \
    reduction(+: approx)
    for (int i = 1; i < n; i++)
        approx += f(a + i*h);
    approx = approx*h;

    return approx;
}

코드에서는 approx를 reduction variable로 사용했습니다.

전체 코드는 마찬가지로 아래에서 확인하실 수 있습니다.

https://github.com/junstar92/parallel_programming_study/blob/master/OpenMP/05_omp_trap4.c

GitHub - junstar92/parallel_programming_study

 

parallel for 디렉티브를 사용하면 한 개의 큰 루프를 병렬화하거나 많은 수의 for 루프를 각 루프 전에 여러 개의 parallel for 디렉티브를 연속으로 추가하여 병렬화 할 수 있습니다. 하지만 parallel for 디렉티브를 사용할 때에는 몇 가지 주의사항이 있습니다. 

• for 루프만을 병렬화함 (while이나 do-while은 병렬화하지 않음)
• 반복 횟수가 정해져 있는 for 루프만을 병렬화함
  예를 들어, 무한 루프 for (; ;)나 루프 안에서 break;가 존재하는 경우에는 병렬화할 수 없음

사실 OpenMP는 정형화된 형태의 for 루프만을 병렬화할 수 있습니다.

• 변수 index는 정수형이거나 포인터형이어야 함(ex, float형은 불가)
• start, end, 그리고 incr는 호환가능한 타입이어야 함
• start, end, incr는 루프를 실행하는 동안 변경되어서는 안됨
• 루프가 실행되는 동안, 변수 index는 for문의 "증가 표현식"에 의해서만 수정될 수 있음

병렬화 가능한 for문의 형태

이러한 제약은 시스템이 루프를 시작하기 전에 반복 횟수를 알도록 해줍니다. 유일한 예외 사항은 루프의 본문 코드에서 exit를 호출하는 경우밖에 없습니다.