Contents
- Pthread, OpenMP에서의 행렬 곱 연산 + 전치 행렬(transpose matrix) 사용
- OpenCV library mat을 사용한 행렬 곱 연산
- CUDA libarary로 구현한 행렬 곱 연산
Matrix Multiplication
이번 포스팅에서는 행렬 곱 연산을 여러 방법으로 구현하고, 그 성능을 비교해볼 예정입니다.
WSL2 - ubuntu 20.04 환경에서 실행했으며, 사용한 라이브러리는 4 종류(Pthreads, OpenMP, OpenCV, CUDA)를 사용하여 행렬 곱 연산을 구현해보았습니다. 행렬 곱셈을 위해서 캐논 알고리즘과 같은 특별한 알고리즘은 사용하지 않습니다.
pthread와 OpenMP에서의 행렬 곱은 3중 for문으로 구현되어 있는데, 캐시 미스에 따라 성능이 달라지기 때문에 2가지 방법을 사용했습니다. 이는 밑에서 코드를 보면서 다시 말씀드리겠습니다.
사용한 컴퓨터 사양은 다음과 같습니다.
CPU : 라이젠 7 5800X 8-core
RAM : 32GB
GPU : RTX 3080
Pthreads Libarary
https://github.com/junstar92/parallel_programming_study/blob/master/pthread/13_pth_mat_mul.c
GitHub - junstar92/parallel_programming_study
Contribute to junstar92/parallel_programming_study development by creating an account on GitHub.
github.com
먼저 pthread로 구현한 행렬 곱 코드입니다. 전체 코드는 위의 링크를 참조해주시기 바랍니다.
첫 번째 방법은 일반적으로 캐시를 신경쓰지 않고, 행렬 곱 공식을 그대로 적용한 3중 for문의 코드입니다.
void* Pth_mat_mul1(void* rank)
{
long my_rank = (long)rank;
int local_m = m / thread_count;
int first_row = my_rank * local_m;
int last_row = first_row + local_m;
#ifdef DEBUG
printf("Thread %ld > local_m = %d\n", my_rank, local_m);
#endif
double temp;
for (int i = first_row; i < last_row; i++) {
for (int j = 0; j < k; j++) {
temp = 0.0;
for (int l = 0; l < n; l++) {
temp += A[i*n + l] * B[l*k + j];
}
C[i*k + j] = temp;
}
}
return NULL;
}다만 행 block distribution을 사용하여, 스레드에 행 단위로 연산을 할당하여 계산을 수행합니다. 만약 총 8개의 스레드로 4096 x 4096의 차원을 가지는 행렬 A와 B를 곱셈 연산할 때, 각 스레드는 결과를 담는 행렬인 C에서 512 x 4096 차원의 연산을 할당받습니다.
하지만 위의 코드를 그대로 사용하면, 16 line에서 캐시 미스가 많이 발생합니다. 캐시는 현재 액세스한 메모리와 가까운 메모리에 액세스할 것이라고 가정하여 캐시 라인에 인접한 메모리 블록을 올려둡니다. 가장 안쪽의 for문에서 A의 경우에는 \(l\)이 증가하면서 바로 옆에 있는 메모리(1열씩 옮기며)에 액세스하지만, B의 경우에는 바로 옆에 있는 인접한 메모리가 아닌 k만큼 떨어진 위치의 메모리, 즉, 행을 한칸씩 내려가면서 메모리에 액세스합니다. 따라서 B 행렬에서 캐시 미스가 계속해서 발생하여 메인 메모리에 액세스하는 시간 때문에 성능이 좋지 못합니다.
위 그림처럼 (4, 4) 위치의 값을 구하기 위해서는 행렬 A의 4행의 요소들과 행렬 B의 4열의 요소들을 곱해야하는데, 행렬 B의 4열 요소들에 접근할 때마다 캐시미스가 발생하게 됩니다.
이 캐시 미스를 줄이기 위해서는 결국 for 반복문을 돌면서 행 단위로 액세스하는 것이 아닌 메모리상 바로 옆에 인접한 열 단위로 액세스해야합니다. 따라서 행렬 A와 행렬 B를 곱할 때, 행렬 B의 행과 열을 바꿔주어 행렬 B의 전치행렬을 사용하기로 결정했습니다.
위와 같이 행렬 B의 행과 열을 바꿔주어 B의 전치 행렬을 사용하게 되면, (4, 4) 위치의 값을 구할 때, 행렬 A의 4행의 요소들과 전치 행렬 B의 4행의 요소들을 곱해주면 됩니다. 그럼 행렬 B에서 행 단위로 메모리에 액세스하던 것이 열 단위로 바로 인접한 메모리에 접근할 수 있게 됩니다. 다만, 전치 행렬 B를 구하기 위해서 행렬 B 크기만큼 메모리가 더 필요하고, 전치 행렬 B를 구하는 추가 작업이 필요합니다.
void Transpose_matrix(double mat[], double mat_t[], int m, int n)
{
for (int i = 0; i < m; i++) {
for (int j = 0; j < n; j++) {
mat_t[j*m + i] = mat[i*n + j];
}
}
}전치 행렬은 위의 코드로 구할 수 있습니다. 그리고 전치 행렬을 사용한 행렬 곱 함수는 다음과 같습니다.
void* Pth_mat_mul2(void* rank)
{
long my_rank = (long)rank;
int local_m = m / thread_count;
int first_row = my_rank * local_m;
int last_row = first_row + local_m;
#ifdef DEBUG
printf("Thread %ld > local_m = %d\n", my_rank, local_m);
#endif
double temp;
for (int i = first_row; i < last_row; i++) {
for (int j = 0; j < k; j++) {
temp = 0.0;
for (int l = 0; l < n; l++) {
temp += A[i*n + l] * BT[j*n + l];
}
C[i*k + j] = temp;
}
}
return NULL;
}16 line에서 A와 BT 모두 for문을 반복하면서 바로 옆에 인접한 메모리에 계속해서 액세스하게 됩니다.
main 함수는 다음과 같습니다.
int main(int argc, char* argv[])
{
Get_args(argc, argv);
pthread_t* thread_handles = (pthread_t*)malloc(thread_count * sizeof(pthread_t));
A = (double*)malloc(m * n * sizeof(double));
B = (double*)malloc(n * k * sizeof(double));
C = (double*)malloc(m * k * sizeof(double));
BT = (double*)malloc(k * n * sizeof(double));
Generate_matrix(A, m, n);
Generate_matrix(B, n, k);
#ifdef DEBUG
Print_matrix(A, m, n, "A");
Print_matrix(B, m, n, "B");
#endif
void* sol_function;
switch (sol) {
case 1:
sol_function = &Pth_mat_mul1;
break;
case 2:
sol_function = &Pth_mat_mul2;
break;
}
double start, finish, avg_elapsed = 0.0;
for (int count = 0; count < NCOUNT; count++) {
GET_TIME(start);
if (sol == 2) {
Transpose_matrix(B, BT, n, k);
}
for (long thread = 0; thread < thread_count; thread++)
pthread_create(&thread_handles[thread], NULL, sol_function, (void*)thread);
for (long thread = 0; thread < thread_count; thread++)
pthread_join(thread_handles[thread], NULL);
GET_TIME(finish);
printf("[%3d] Elapsed time = %.6f seconds\n", count+1, finish-start);
avg_elapsed += (finish - start) / NCOUNT;
}
#ifdef DEBUG
Print_matrix(C, m, k, "The product is");
#endif
printf("Average elapsed time = %.6f seconds\n", avg_elapsed);
free(A);
free(B);
free(C);
free(BT);
free(thread_handles);
return 0;
}프로그램을 실행할 때, 스레드의 개수(thread_count), 행렬 차원(m, n, k), 그리고 행렬 곱 계산 방법을 입력으로 전달받습니다. 1이 캐시 미스를 고려하지 않은 일반적인 방법이며, 2가 행렬 B의 전치 행렬을 구해 캐시 미스를 줄이는 방법입니다. 성능 측정을 위해서 동일한 연산을 100번 수행하여 평균 시간을 출력합니다. 매 iteration마다 수행 시간을 출력하기는 합니다. 그리고 2번 방법에 대해서 수행 시간에 행렬 B를 전치하는 것도 포함하였습니다.
컴파일은 다음 커맨드로 하고,
$ gcc -Wall -o 13_pth_mat_mul 13_pth_mat_mul.c -pthread16개의 스레드와 4096x4096 행렬 A, B의 곱을 1번 방법으로 실행한 결과는 다음과 같습니다.
2번 방법으로 실행한 결과입니다.
전치 행렬을 구하는 작업이 추가됬음에도 거의 2배 정도 성능이 향상된 것을 확인할 수 있습니다.
OpenMP Libarary
https://github.com/junstar92/parallel_programming_study/blob/master/OpenMP/14_omp_mat_mul.c
GitHub - junstar92/parallel_programming_study
Contribute to junstar92/parallel_programming_study development by creating an account on GitHub.
github.com
이번에는 OpenMP로 구현한 코드입니다. 전체 코드는 위 링크를 참조해주세요.
OpenMP는 pthread와 구현 자체는 비슷합니다만, 당연히 OpenMP 라이브러리를 사용하여 멀티 스레드 환경으로 행렬 곱을 수행합니다. 또한, pthread와 동일하게 1번 방법과 2번 방법을 모두 적용하여 비교할 수 있습니다. 한 가지 차이점은 행렬 B의 전치 행렬을 멀티 스레드 환경으로 구합니다. 따라서 전치 행렬을 구하는 작업은 pthread보다 조금 더 빠를 것입니다.
먼저 캐시 미스를 고려하지 않은 3중 for문으로 행렬의 곱을 구하는 함수입니다.
void Omp_mat_mul1(double A[], double B[], double C[], int m, int n, int k, int thread_count)
{
double start, finish, temp, avg_elapsed = 0.0;
for (int count = 0; count < NCOUNT; count++) {
start = omp_get_wtime();
#pragma omp parallel for num_threads(thread_count) \
default(none) private(temp) shared(A, B, C, m, n, k)
for (int i = 0; i < m; i++) {
for (int j = 0; j < k; j++) {
temp = 0.0;
for (int l = 0; l < n; l++) {
temp += A[i*n + l] * B[l*k + j];
}
C[i*k + j] = temp;
}
}
finish = omp_get_wtime();
printf("[%3d] Elapsed time = %f seconds\n", count+1, finish - start);
avg_elapsed += (finish-start) / NCOUNT;
}
printf("Average elapsed time : %.6f seconds\n", avg_elapsed);
}for문에 parallel for 디렉티브가 추가되는 것외에는 pthread와 크게 달라지는 것은 없습니다.
다음은 캐시 미스를 고려하여 B의 전치 행렬을 구해 행렬 곱 연산을 하는 함수입니다.
void Omp_mat_mul2(double A[], double B[], double C[], double BT[], int m, int n, int k, int thread_count)
{
double start, finish, temp, avg_elapsed = 0.0;
for (int count = 0; count < NCOUNT; count++) {
start = omp_get_wtime();
Transpose_matrix(B, BT, n, k, thread_count);
#pragma omp parallel for num_threads(thread_count) \
default(none) private(temp) shared(A, BT, C, m, n, k)
for (int i = 0; i < m; i++) {
for (int j = 0; j < k; j++) {
temp = 0.0;
for (int l = 0; l < n; l++) {
temp += A[i*n + l] * BT[j*n + l];
}
C[i*k + j] = temp;
}
}
finish = omp_get_wtime();
printf("[%3d] Elapsed time = %f seconds\n", count+1, finish - start);
avg_elapsed += (finish-start) / NCOUNT;
}
printf("Average elapsed time : %.6f seconds\n", avg_elapsed);
}pthread와 마찬가지로 line 13에서 캐시 미스를 최소화하기 위해 매 반복마다 인접한 메모리에 액세스할 수 있도록 하였습니다.
gcc -Wall -fopenmp -o 14_omp_mat_mul 14_omp_mat_mul.c위 커맨드로 컴파일 한 후에, 1번 방법과 2번 방법으로 실행한 결과입니다.
pthread와 OpenMP에서의 성능 차이는 거의 없었습니다.
OpenCV Libarary
https://github.com/junstar92/parallel_programming_study/blob/master/ocv_mat_mul.c
GitHub - junstar92/parallel_programming_study
Contribute to junstar92/parallel_programming_study development by creating an account on GitHub.
github.com
OpenCV에는 mat 타입을 제공합니다. 그리고 두 mat 타입을 곱하면 알아서 행렬 곱 연산을 수행해주는데, 이 연산의 수행 시간을 측정하였습니다. 전체 코드는 위 링크를 참조해주세요.
측정을 위한 main 함수는 다음과 같습니다.
#include <stdio.h>
#include <stdlib.h>
#include <opencv2/highgui.hpp>
#include <sys/time.h>
int main(int argc, char* argv[])
{
int m, n, k;
Get_args(argc, argv, &m, &n, &k);
double *A, *B, *C;
A = (double*)malloc(m * n * sizeof(double));
B = (double*)malloc(n * k * sizeof(double));
C = (double*)malloc(m * k * sizeof(double));
Generate_matrix(A, m, n);
Generate_matrix(B, n, k);
#ifdef DEBUG
Print_matrix(A, m, n, "A");
Print_matrix(B, m, n, "B");
#endif
double start, finish, avg_elapsed = 0.0;
for (int count = 0; count < NCOUNT; count++) {
GET_TIME(start);
cv::Mat cvA(m, n, CV_64FC1, A);
cv::Mat cvB(n, k, CV_64FC1, B);
cv::Mat cvC = cvA * cvB;
//cv::gemm(cvA, cvB, 1.0, NULL, 0, cvC);
//C = reinterpret_cast<double*>(cvC.data);
GET_TIME(finish);
printf("[%3d] Elapsed time = %.6f seconds\n", count+1, finish-start);
avg_elapsed += (finish - start) / NCOUNT;
}
#ifdef DEBUG
printf("The product is\n");
cv::print(C);
printf("\n\n");
#endif
printf("Average elapsed time = %.6f seconds\n", avg_elapsed);
free(A);
free(B);
free(C);
return 0;
}최대한 pthread와 OpenMP 때와 동일한 환경으로 비교하기 위해서 double 타입으로 메모리를 할당하고, 행렬 A, B값을 초기화한 다음에 곱 연산을 수행하기 전에 cv::Mat 타입으로 각 행렬을 선언해줍니다. 핸디캡처럼 cv::Mat 타입을 생성하는 시간까지 수행 시간에 포함하였으며(line 27-28), 곱 연산은 line 29에서 이루어집니다.
pthread와 OpenMP와는 달리 행렬의 크기(m, n, k)만 입력으로 받습니다.
사용하는 OpenCV 버전에 따라서 컴파일 커맨드가 조금 달라질 수는 있겠지만, 저의 경우에는 OpenCV 4.2.0을 사용하고 있고 아래의 커맨드로 컴파일을 수행하였습니다. 'pkg-config opencv4 --libs --cflags'는 opencv 라이브러리를 사용하고 링크하기 위한 플래그로 대체됩니다.
g++ -Wall -o ocv_mat_mul ocv_mat_mul.c $(pkg-config opencv4 --libs --cflags)
실행 결과는 다음과 같습니다.
pthread와 OpenMP와 비교해서 캐시 미스를 신경쓰지 않은 1번 방법보다는 2배가량 좋은 성능을 보여주고 있지만, 캐시 미스를 고려한 2번 방법보다는 아주 약간 성능이 떨어지는 모습을 보여줍니다. OpenCV 내부에서 행렬 곱셈이 어떻게 구현되어 있는지는 알아봐야겠지만, 단순히 위의 방법으로는 (2번 방법을 사용한)pthread와 OpenMP보다 아주 약간 성능이 떨어지는 것을 보여주고 있습니다.
CUDA Library
https://github.com/junstar92/parallel_programming_study/blob/master/cuda_mat_mul.cu
GitHub - junstar92/parallel_programming_study
Contribute to junstar92/parallel_programming_study development by creating an account on GitHub.
github.com
GPU를 사용하여 행렬 곱을 병렬로 연산하도록 구현하였습니다. 전체 코드는 위를 참조해주시기 바랍니다. 단순히 행렬 곱셈 성능을 비교하기 위해서 CUDA 라이브러리를 사용하여 구현해봤는데, CUDA와 관련한 자세한 내용은 이번 글에서는 생략하려고 합니다. GPU와 CUDA에 관해서는 나중에 따로 스터디하여 포스팅해보도록 하겠습니다.
main 함수와 GPU device에서 수행되는 커널 함수만 가져와봤습니다.
#include <stdio.h>
#include <stdlib.h>
#include <sys/time.h>
#include <cuda_runtime.h>
__global__ void cuda_mat_mul(double *A, double *B, double *C, int m, int n, int k);
int main(int argc, char* argv[])
{
int m, n, k;
Get_args(argc, argv, &m, &n, &k);
double *A, *B, *C;
A = (double*)malloc(m * n * sizeof(double));
B = (double*)malloc(n * k * sizeof(double));
C = (double*)malloc(m * k * sizeof(double));
Generate_matrix(A, m, n);
Generate_matrix(B, n, k);
#ifdef DEBUG
Print_matrix(A, m, n, "A");
Print_matrix(B, n, k, "B");
#endif
// Allocate the device input matrixs for A, B, C;
double* d_A, *d_B, *d_C;
CUDA_CHECK(cudaMalloc((void**)&d_A, m * n * sizeof(double)));
CUDA_CHECK(cudaMalloc((void**)&d_B, n * k * sizeof(double)));
CUDA_CHECK(cudaMalloc((void**)&d_C, m * k * sizeof(double)));
dim3 dimGrid(ceil(n / 32.0), ceil(m / 32.0));
dim3 dimBlock(32, 32); // 256 threads
double start, finish, avg_elapsed = 0.0;
// Launch the Matrix Multiplication CUDA Kernel
CUDA_CHECK(cudaSetDevice(0));
for (int count = 0; count < NCOUNT; count++) {
GET_TIME(start);
// Copy the host matrixs A and B in host memory to the device matrixs in device memory
CUDA_CHECK(cudaMemcpy(d_A, A, m * n * sizeof(double), cudaMemcpyHostToDevice));
CUDA_CHECK(cudaMemcpy(d_B, B, n * k * sizeof(double), cudaMemcpyHostToDevice));
cuda_mat_mul<<<dimGrid, dimBlock>>>(d_A, d_B, d_C, m, n, k);
CUDA_CHECK(cudaDeviceSynchronize());
CUDA_CHECK(cudaGetLastError());
// Copy the device result matrix in device memory to the host result matrix in host memory
CUDA_CHECK(cudaMemcpy(C, d_C, m * k * sizeof(double), cudaMemcpyDeviceToHost));
GET_TIME(finish);
printf("[%3d] Elapsed time = %.6f seconds\n", count+1, finish-start);
avg_elapsed += (finish - start) / NCOUNT;
}
#ifdef DEBUG
Print_matrix(C, m, k, "The product is");
#endif
printf("Average elapsed time = %.6f seconds\n", avg_elapsed);
// Free device global memory
CUDA_CHECK(cudaFree(d_A));
CUDA_CHECK(cudaFree(d_B));
CUDA_CHECK(cudaFree(d_C));
free(A);
free(B);
free(C);
return 0;
}
__global__ void cuda_mat_mul(double *A, double *B, double *C, int m, int n, int k)
{
int ROW = blockIdx.x * blockDim.x + threadIdx.x;
int COL = blockIdx.y * blockDim.y + threadIdx.y;
if (ROW < m && COL < k) {
double value = 0.0;
for (int i = 0; i < k; i++) {
value += A[ROW * n + i] * B[i * k + COL];
}
C[ROW * k + COL] = value;
}
}GPU에서 연산을 하기 위해서는 CPU 메모리에 있는 행렬 값들을 GPU 메모리로 복사하는 과정이 필요한데, 이 과정을 수행 시간에 포함하였습니다. 그리고 곱 연산 결과가 GPU 메모리에 있기 때문에, 이 결과를 CPU 메모리로 복사하는 과정 또한 수행 시간에 포함했습니다. 최대한 비슷한 환경을 맞춘다고는 했는데, 아직 CUDA에 대해서 많이 알지 못하기 때문에 부족한 부분이 있을 수 있습니다... !
다음의 커맨드로 컴파일을 할 수 있습니다.
nvcc -o cuda_mat_mul cuda_mat_mul.cu $(pkg-config cuda-XX.X --libs --cflags)
# or nvcc -o cuda_mat_mul cuda_mat_mul.cu -lcudaOpenCV와 마찬가지로 pkg-config를 통해서 CUDA 라이브러리의 컴파일/링킹 옵션을 가져왔습니다. cuda-XX.X는 본인에게 맞는 버전을 지정해주면 됩니다. 저는 CUDA 11.1 버전을 사용하고 있습니다.
(pkg-config 대신 '-lcuda' 를 사용하도 가능합니다.)
실행 결과입니다.
역시 GPU의 성능은 최고입니다...
CUDA에 대한 지식이 깊지 않아서 최적화도 되지 않은 아주 기본적인 코드이지만 다른 라이브러리에서 아무리 빨라도 14초~15초가 걸리던 연산을 단 0.6초 정도만에 완료하고 있습니다.
pthread와 OpenMP에서 특별한 알고리즘을 사용하지 않고 행렬 곱 연산을 구현해봤는데, OpenCV와 큰 차이가 없다는 점에서 조금 놀랐지만.. 16개의 스레드를 사용했을 때, CPU 사용량이 거의 100%를 찍고 있습니다.
반면에 OpenCV의 경우에는 약 10~12%의 CPU 사용량을 보여주고 있어서, OpenCV의 최적화가 어떻게 되어있는지 궁금하네요.. 대강 구글링했을 때 찾지는 못했는데, 코드를 뜯어보던가 해야겠습니다... ㅠ
(하지만 결론은 'GPU 최고'입니다)
아래 이미지는 각 라이브러리 케이스 별로 비교하기 쉽도록 1회만 수행한 결과입니다.
엄격하게 동일한 환경을 구성하고 구현한 것이 아니기 때문에 수행 결과는 그냥 참고용으로만 보시는게 좋습니다.. !
이상입니다 !
'프로그래밍 > 병렬프로그래밍' 카테고리의 다른 글
| 행렬 곱 연산 비교 2 (OpenBLAS, Intel MKL, cuBLAS) (0) | 2021.11.27 |
|---|---|
| [OpenMP] 행렬-벡터 곱 연산 (0) | 2021.11.25 |
| [OpenMP] Critical Section (0) | 2021.11.25 |
| [OpenMP] schedule clause (0) | 2021.11.24 |
| [OpenMP] for 디렉티브 (0) | 2021.11.23 |
댓글