[MPI] 행렬 - 벡터 곱 연산 + 성능 평가

References

• An Introduction to Parallel Programming

Contents

• 행렬 - 벡터 곱 연산
• MPI_Allgather
• MPI_Wtime, MPI_Barrier

---

이번에는 행렬-벡터 곱 연산을 MPI를 사용하여 병렬화해보도록 하겠습니다.

 

\(\text{A}\)가 m x n 행렬이고, \(\textbf{x}\)가 n개의 컴포넌트를 갖고 있는 벡터라면, \(\textbf{y} = \text{A}\textbf{x}\)는 m개의 컴포넌트를 갖는 벡터가 됩니다. 그리고 y의 i번째 컴포넌트는 행렬 A의 i번째 행과 x의 dot product로 계산할 수 있습니다.

 

이 연산을 시리얼로 의사 코드를 작성하면, 다음과 같습니다.

위 코드는 일반적인 C코드입니다. 위에서 행렬을 표현하기 위해서 2차원배열인 A를 사용했는데, 아래에서 실제 구현할 때에는 1차원배열을 2차원배열처럼 사용할 것입니다. 잠시 뒤에 살펴볼 예제에서는 다음과 같은 2차원배열을 

아래의 1차원배열로 저장합니다.

따라서, 0부터 행과 열을 카운트하는데, 2차원배열의 2행 1열에 저장되어 있는 요소(9)은 \(2 \times 4 + 1 = 9\)로 구할 수 있습니다. 좀 더 일반적으로 표현하면 배열이 n개의 열을 갖는다면 i번째 행과 j번째 열에 저장되어 있는 요소는 일차원 배열의 \(i \times n + j\)의 위치에 저장됩니다. 이 1차원 배열을 2차원 배열처럼 사용하는 방식으로 행렬 - 벡터 곱 함수를 작성하면 아래처럼 작성할 수 있습니다.

void Mat_vec_mul(
	double A[]  /* in  */,
    double B[]  /* in  */,
    double C[]  /* out */,
    int    m    /* in  */,
    int    n    /* in  */)
{
	for (int i = 0; i < m; i++) {
    	y[i] = 0.0;
        for (int j = 0; j < n; j++)
        	y[i] += A[i*n + j] * x[j];
    }
}

 

이제 이 함수를 어떻게 병렬화하는지 알아보겠습니다. 개별 태스크는 A의 요소와 x의 컴포넌트의 곱과 이 결과를 y의 컴포넌트에 더하는 것입니다. 즉, 아래 문장의 실행이

\[y[i] += A[i*n + j] * x[j];\]

하나의 태스크가 됩니다.

따라서 만약 \(y[i]\)를 프로세스 q에 할당한다면, 이는 A의 i행을 프로세스 q에 할당하는 것으로 편리하게 태스크를 각 프로세스에 할당할 수 있습니다.

즉, block distribution, cyclic distribution, 또는 block-cyclic distribution을 사용하여 A를 행(row)으로 할당할 수 있습니다. MPI에서는 block distribution을 사용하는 것이 가장 쉽기 때문에 우리는 A의 행을 block distribution을 사용하여 파티셔닝할 것이고 보통의 경우, m행을 comm_sz로 나누어서 m/comm_sz 행씩 각 프로세스에 할당합니다.

 

행으로 A를 분산하고 \(y[i]\)의 연산이 A의 요소를 포함하기 때문에 y 또한 block으로 분산되어야 합니다. 즉, A의 i번째 행이 프로세스 q로 할당되면, y의 i번째 컴포넌트 또한 프로세스 q에 할당되어야 합니다.

 

\(y[i]\)의 연산은 A의 i번째 행과 x의 모든 컴포넌트를 포함하기 때문에 통신의 양은 x 전체를 각 프로세스에 할당하는 것만으로 통신을 최소화할 수 있습니다. 따라서 x가 block distribution으로 분산될 때, 각 프로세스가 x의 전체 컴포넌트에 어떻게 액세스하는 지만 고민하면 됩니다.

이전 게시글(데이터 분산)에서 설명한 트리 구조의 브로드캐스트를 사용하는 MPI_Bcast와 MPI_Gather를 사용하여 접근할 수 있습니다. 하지만, 이런 경우에는 버터플라이(butterfly) 구조의 트리를 사용하는 MPI_Allgather 함수를 사용하는 것이 좋습니다.

int MPI_Allgather(
	void*        send_buf_p /* in  */,
    int          send_count /* in  */,
    MPI_Datatype send_type  /* in  */,
    void*        recv_buf_p /* out */,
    int          recv_count /* in  */,
    MPI_Datatype recv_type  /* in  */,
    MPI_Comm     comm       /* in  */);

이 함수는 각 프로세스의 send_buf_p에 저장된 값들을 각 프로세스의 recv_buf_p가 가리키는 메모리에 저장합니다. 보통 recv_count는 각 프로세스에서 전달받은 데이터의 크기이고, recv_count는 send_count와 같습니다.

 

이제 위의 시리얼 행렬-벡터 곱 함수처럼 병렬 행렬-벡터 곱 함수를 작성하면, 다음과 같이 작성할 수 있습니다.

void Mat_vec_mul(
    double      local_A[]   /* in  */,
    double      local_x[]   /* in  */,
    double      local_y[]   /* out */,
    int         local_m     /* in  */,
    int         n           /* in  */,
    int         local_n     /* in  */,
    MPI_Comm    comm        /* in  */)
{
    double* x;
    int local_ok = 1;

    x = (double*)malloc(n * sizeof(double));
    if (x == NULL)
        local_ok = 0;
    Check_for_error(local_ok, "Mat_vec_mul",
        "Can't allocate temporary vector", comm);
    
    MPI_Allgather(local_x, local_n, MPI_DOUBLE,
        x, local_n, MPI_DOUBLE, comm);
    
    for (int local_i = 0; local_i < local_m; local_i++) {
        local_y[local_i] = 0.0;
        for (int j = 0; j < n; j++)
            local_y[local_i] += local_A[local_i * n + j] * x[j];
    }
    free(x);
}

실행해보면 위의 출력 결과와 같이 나오게 됩니다.

전체 코드는 아래 링크에서 확인 가능합니다.

https://github.com/junstar92/parallel_programming_study/blob/master/mpi/08_mpi_mat_vec_mul.c

GitHub - junstar92/parallel_programming_study

---

 

성능 측정 (MPI_Wtime)

보통 성능은 연산에 걸리는 시간을 측정하여 평가합니다. 그리고, 그 결과를 출력할 때걸리는 시간이나 매트릭스에 데이터를 입력하는 시간은 관심이 없기 때문에 수행 시간은 프로그램이 시작해서 끝날 때까지가 아닌 연산에 걸리는 시간만 측정합니다. 따라서, 실제 행렬-벡터 곱셈 연산의 시작부터 끝까지 걸리는 시간을 재는 함수를 소스 코드에 추가해야 합니다.

 

MPI는 시간을 측정하기 위한 함수로 MPI_Wtime을 제공하며, 이 함수는 과거 어느 시점부터의 걸린 시간(초단위)을 리턴합니다.

double MPI_Wtime(void);

따라서 다음과 같이 MPI 코드의 특정 부분의 시간을 잴 수 있습니다.

double start, finish;
...
start = MPI_Wtime();
/* Code to be timed */
...
finish = MPI_Wtime();
printf("Proc %d > Elapsed time = %f seconds\n", my_rank, finish-start);

 

행렬-벡터 곱 병렬 연산과 시리얼 연산을 시간 측정을 통해서 비교해볼텐데, 시리얼 연산에서는 MPI 라이브러리를 링크할 필요없이, gettimeofday라는 POSIX 라이브러리 함수를 호출하여 시간을 측정할 수 있습니다. 이 함수는 과거 어느 시점부터 걸린 시간(마이크로초단위)을 리턴합니다. 

이는 아래 정의한 매크로를 사용하여 다음과 같이 측정할 수 있습니다. 매크로를 실행하고 나면 now에는 초 단위의 시간이 저장됩니다.

#include <sys/time.h>

/* The argument now should be a double (not a pointer to a double) */
#define GET_TIME(now) \
{ \
    struct timeval tv; \
    gettimeofday(&tv, NULL); \
    now = tv.tv_sec + tv.tv_usec/1000000.0; \
}

...

double start, finish;
...
GET_TIME(start);
/* Code to be timed */
...
GET_TIME(finish);
printf("Elapsed time = %f seconds\n", finish-start);

time.h 헤더에서는 clock 함수를 제공하는데 이는 CPU 시간을 리턴합니다. CPU 시간에는 유저 코드, 라이브러리 함수, OS 코드에서 소비된 모든 시간을 포함하고, idle 시간(sleep에 빠져있었던 시간)은 포함하지 않습니다. 병렬 실행에서는 이 idle시간도 중요한 부분이므로 병렬 프로그램에서 clock 함수는 적절하지 않습니다. 반면에 gettimeofday나 MPI_Wtime은 절대 시간을 리턴합니다.

 

병렬 프로그램에서 살펴볼 것이 한 가지 더 있는데, 아시다시피 병렬 프로그램은 comm_sz개의 프로세스가 동시에 수행되지만 우리가 원하는 것은 하나의 시간입니다. 이상적으로는 모든 프로세스가 행렬-벡터 곱셈 연산을 하여 마지막 프로세스가 끝났을 때 걸리는 시간을 리포트하면 됩니다. 다시 말하면 병렬 프로그램에서 수행 시간은 "가장 느린" 프로세스가 끝날 때까지 걸린 시간이 됩니다. 모든 프로세스가 동시에 시작하는 것을 보장하지 못하면 이 시간을 정확하게 구할 수가 없습니다. 이를 위해 MPI 집합 통신 함수 MPI_Barrier를 제공하는데, 이 함수는 모든 프로세스가 이 함수를 호출할 때까지 어떠한 프로세스도 리턴하지 않도록 보장해줍니다. 

int MPI_Barrier(MPI_Comm comm /* in */);

다음의 코드의 line 21 ~ 30은 병렬 프로그램에서 MPI 코드 블록의 시간을 재는 코드이며, 각 프로세스에서 걸린 시간들을 하나의 시간으로 리포트해줍니다.

int main(void)
{
    double *local_A, *local_x, *local_y;
    int m, local_m, n, local_n;
    int my_rank, comm_sz;
    MPI_Comm comm;

    MPI_Init(NULL, NULL);
    comm = MPI_COMM_WORLD;
    MPI_Comm_size(comm, &comm_sz);
    MPI_Comm_rank(comm, &my_rank);

    Get_dims(&m, &local_m, &n, &local_n, my_rank, comm_sz, comm);
    Allocate_arrays(&local_A, &local_x, &local_y, local_m, n, local_n, comm);

    Get_matrix(local_A, m, local_m, n, "A", my_rank, comm);
    //Print_matrix(local_A, m, local_m, n, "A", my_rank, comm);
    Get_vector(local_x, n, local_n, "x", my_rank, comm);
    //Print_vector(local_x, n, local_n, "x", my_rank, comm);

    double local_start, local_finish, local_elapsed, elapsed;
    local_start = MPI_Wtime();
    Mat_vec_mul(local_A, local_x, local_y, local_m, n, local_n, comm);
    local_finish = MPI_Wtime();
    local_elapsed = local_finish - local_start;
    MPI_Reduce(&local_elapsed, &elapsed, 1, MPI_DOUBLE, MPI_MAX, 0, comm);
    //Print_vector(local_y, m, local_m, "y", my_rank, comm);

    if (my_rank == 0)
        printf("Elapsed time = %f seconds\n", elapsed);

    free(local_A);
    free(local_x);
    free(local_y);

    MPI_Finalize();

    return 0;
}

 

그리고 행렬-벡터 시리얼 곱셈은 다음의 코드로 실행합니다.

int main(void)
{
    double *A, *x, *y;
    int m, n;

    Get_dims(&m, &n);
    Allocate_arrays(&A, &x, &y, m, n);

    Get_matrix(A, m, n, "A");
    //Print_matrix(A, m, n, "A");
    Get_vector(x, n, "x");
    //Print_vector(x, n, "x");

    double start, finish;
    GET_TIME(start);
    Mat_vec_mul(A, x, y, m, n);
    GET_TIME(finish);
    Print_vector(y, m, "y");

    printf("Elapsed time = %f seconds\n", finish - start);

    free(A);
    free(x);
    free(y);

    return 0;
}

 

각 시리얼과 병렬 코드의 전체 내용은 아래 링크에서 확인 가능합니다 !

https://github.com/junstar92/parallel_programming_study/blob/master/mpi/09_serial_mat_vec_mul_time.c

https://github.com/junstar92/parallel_programming_study/blob/master/mpi/10_mpi_mat_vec_mul_time.c

 

각 코드를 실행하면, 다음과 같이 수행 시간이 출력됩니다.

시리얼 곱셈 연산 결과

병렬 곱셈 연산 결과

16x16의 행렬과 벡터의 곱셈은 크기가 작아서 시리얼 연산의 속도가 더 빠르게 나오는 것을 확인할 수 있습니다.

 

다음은 행렬의 크기와 프로세스의 개수에 따라 측정한 결과입니다. 행렬의 크기는 m과 n을 같은 값으로 두고 평가하였습니다.

행렬-벡터 병렬 곱셈 연산 결과 (단위 : ms)

행렬-벡터 시리얼 곱셈 연산 결과 (단위 : ms)

병렬 곱셈 연산 결과에서 process가 1인 경우에는 싱글 코어에서 동작하는 시리얼 프로그램의 실행 시간과 거의 동일하게 나옵니다. process의 수를 고정시키고 매트릭스의 차수를 2배씩 증가시키면 실행 시간은 4배 정도 증가하는 것을 볼 수 있습니다. 

매트릭스 차수를 고정시키고 process의 수를 증가시키면 실행 시간은 감소하는 모습을 보여주고 있습니다. 프로세스의 수가 더 많으면 어떻게 여기서 알 수는 없지만, 대체로 증가하다가 어느 순간부터 전체 실행 시간은 변화가 없을 것으로 예상합니다. 또한, 프로세스가 4개까지는 2배 커질 때 걸리는 시간도 2배로 감소하다가 그 이후부터는 2배 이하로 증가하는 모습을 보여주고 있습니다.