[Multi-Core & GPU Programming] HW2 - Matrix Verification Challenge
"Matrix Verification Challenge"
[Objective]
GEMM 방식과, GEMV 방식을 사용한 Freivalds'Algorithm을 사용하여 "A x B == C"를 검증
[Theory]
▣ GEMM(General Matrix-Matrix Multiplication)
· 행렬과 행렬의 곱셈으로 N x N 행렬 2개를 곱하려면 가로 열과 세로 열을 모두 곱하고 더해야 하므로, 시간 복잡도가 O(N3)
→ 데이터 크기가 커지면 연산 시간이 기하급수적으로 증가
▣ GEMV(General Matrix-Vector Multiplication)
· 행렬과 Vector의 곱셈으로 행렬의 각 행과 Vector를 곱하기만 하면 되므로 시간 복잡도는 O(N2)
→ GEMM에 비해 훨씬 가볍고 빠름
▣ Freivalds' Algorithm
· 계산해 놓은 행렬 곱셈 결과(C)가 정말 정답(A x B)이 맞는지 초고속으로 검증하는 Algorithm
· A x B를 다시 계산하려면 O(N3) GEMM → 임의의 1차원 Vector를 하나를 만들어서 A x (B x v) = C x v가 성립하는지 확인
▣ Deterministic
· 동일한 입력이 주어지면 100% 동일한 결과와 정답을 보장하는 방식
→ Parallel GEMM 방식처럼 행렬 전체를 직접 다 곱해서 한 치의 오차도 없이 비교
△ 정확도↑
▼ 속도↓
▣ Probabilistic
· 내부적으로 난수를 사용하여 연산 과정을 대폭 생략하는 대신, 아주 미세한 확률로 틀린 답을 낼 가능성을 감수하는 방식
△ 속도↑
▼ 정확도↓
[Implementation]
▣ Code Description
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 | #include <algorithm> #include <cassert> #include <cmath> #include <cstdlib> #include <numeric> #include <thread> // Declare Header To Use Thread #include <random> // Declare Header To Use Random Number /** * @brief Initializes a vector with random double values. * * This function fills the given array * * @param a Pointer to the array to be initialized. * @param N The number of elements in the array. */ inline void init_vec(double *a, int N) { const int nthreads = 32; // # Of Threads int chunk_size = N / nthreads; // Chunk Size std::thread threads[nthreads]; // Array Of Threads for (int t = 0; t < nthreads; t++) { int start = chunk_size * t; // Start Index int end = (t == nthreads - 1) ? N : start + chunk_size; // End Index(Allocate Remainder To Last Thread) threads[t] = std::thread([=]() { std::random_device rand; // To Get Random Seed std::mt19937 gen(rand() + t); // Create Random Number std::uniform_real_distribution<double> dist(0.0, 1.0); // Create Double Value for (int i = start; i < end; i++) { a[i] = dist(gen); // Get One Random Number } }); } for (int t = 0; t < nthreads; t++) { threads[t].join(); // Threads Join } } /** * @brief Performs a matrix-vector multiplication. * * This function computes the product of a matrix 'a' and a vector 'b', storing * the result in vector 'c'. * * @param a Pointer to the first element of the matrix 'a' (assumed to be in * row-major order). * @param b Pointer to the first element of the vector 'b'. * @param c Pointer to the first element of the result vector 'c'. * @param N The dimension of the matrix and vectors (assuming a square matrix * and compatible vector sizes). */ inline void gemv(double *a, double *b, double *c, int N) { std::fill(c, c + N, 0); // Initialize Matrix C To 0 const int nthreads = 32; // # Of Threads int chunk_size = N / nthreads; // Chunk Size std::thread threads[nthreads]; // Array Of Threads for (int t = 0; t < nthreads; t++) { int start = chunk_size * t; // Start Index int end = (t == nthreads - 1) ? N : start + chunk_size; // End Index(Allocate Remainder To Last Thread) threads[t] = std::thread([=]() { for (int i = start; i < end; i++) { for (int j = 0; j < N; j++) { c[i] += a[i * N + j] * b[j]; // C[i] += A[i][j] x B[j] } } }); } for (int t = 0; t < nthreads; t++) { threads[t].join(); // Thread Join } } /** * @brief Performs matrix multiplication of two NxN matrices. * * This function computes the product of two square matrices `a` and `b`, * and stores the result in matrix `c`. All matrices are represented as * 1-dimensional arrays in row-major order. * * @param a Pointer to the first input matrix (NxN). * @param b Pointer to the second input matrix (NxN). * @param c Pointer to the output matrix (NxN) where the result will be stored. * @param N The dimension of the matrices (number of rows and columns). */ inline void gemm(double *a, double *b, double *c, int N) { std::fill(c, c + N * N, 0); // Initialize Matrix C To 0 const int nthreads = 32; // # Of Threads int chunk_size = N / nthreads; // Chunk Size std::thread threads[nthreads]; // Array Of Threads for (int t = 0; t < nthreads; t++) { int start = chunk_size * t; // Start Index int end = (t == nthreads - 1) ? N : start + chunk_size; // Allocate To Last Thread To Prevent Remainder threads[t] = std::thread([=]() { for (int i = start; i < end; i++) { for (int k = 0; k < N; k++) { for (int j = 0; j < N; j++) { c[i * N + j] += a[i * N + k] * b[k * N + j]; // C[i][j] += A[i][k] * B[k][j] } } } }); } for (int t = 0; t < nthreads; t++) { threads[t].join(); // Threads Join } } | cs |
· std::thread([=]( ) { }) : "="으로 캡처하면 Thread가 생성되는 순간의 값이 복사되어 각 Thread는 자기 start/end/t 값을 따로 가짐
> gemm
① 32개의 Thread 생성
② 각 Thread는 자신이 담당할 결과 행렬 C의 행 구간을 할당받음
③ Thread들은 행렬 A, B의 데이터를 읽어와서 곱셈 연산(i → j → k)을 수행하여 행렬 C 생성
④ 완벽한 정답 행렬 C = A x B 생성
** 순서(i → j → k)는 메모리 Access Pattern을 Row-Wise로 변경하여 Cache의 효율성 증가
> init_vec
① 32개의 Thread가 길이 N짜리 빈 Vector v를 N/32개씩 나누어 맡고, 각자 독립적인 난수를 생성하여 무작위 실수로 채움
· gen(rand( ) + t) : "rand( ) + t"는 Thread마다 다른 Seed를 사용해서 서로 다른 난수 패턴이 되도록 함
→ 단순히 "rand( )"인 경우 같은 Seed를 사용하게 되어 같은 난수 패턴 출력 가능성 존재
· uniform_real_distribution : 균등 분포
· dist : "gen"이 만든 난수를 이용해서 0~1 범위의 double 하나를 뽑기
> gemv
① 32개의 Thread를 생성하여 원본 행렬 B와 "init_vec"에서 만든 Random Vector v를 곱하여 임시 Vector 생성
② 행렬 A와 방금 만든 임시 Vector를 곱함
③ "gemm"에서 만든 행렬 C와 Random Vector 를 곱함
④ Thread 연산이 모두 끝난 후, 단일 Thread가 A x (B x v)와 C x v에 들어있는 N개의 숫자를 1번부터 끝까지 1:1로 비교 → O(N)
▣ Result
[Discussion]
· Machine마다 최적 Thread수가 다른데 Thread 수를 32로 고정
→ 실제 HW보다 많은 Thread 생성 : 성능 저하 / 더 많은 Core 보유 : 충분히 활용 X
∴ std::thread::hardware_concurrency( )를 사용하여 Thread 수를 유동적으로 정해서 최적화 가능
· int_vec, gemv, gemm 모두 함수가 호출될 때마다 32개의 Thread를 만들고, 끝나면 join
→ 함수 호출마다 Thread 생성 시 계산보다 Thread 관리 비용 비중↑
∴ Thread Pool(Thread를 한 번만 생성하고 작업만 재사용) / Static Worker(최소한 반복 호출에서 Thread 재활용)로 최적화
· chunk_size의 나머지를 마지막 Thread에 집중
→ 전체 실행 시간은 가장 늦게 끝나는 Thread에 의해 결정되어 마지막 Thread를 기준으로 실행 시간 측정
∴ start = N * t / nthreads, end = N * (t + 1) / nthreads와 같이 균등하게 분배하는 것이 더 효율적
· gemv와 gemm 모두 결과 Vector 전체를 먼저 0으로 초기화한 후 각 Thread에서 결과를 누적
→ 전체를 초기화하므로 시간 지연 발생
∴ 지역 변수를 활용하여 각 Thread마다 초기화한 후 누적 한 다음, 결과값을 한 번만 C 행렬에 Write
→ 메모리 Write 횟수 감소 + Cache 효율 증가 + 실행 시간 감소
[Reference]
· HW2 - Matrix Verification Challenge (MGP) - Yongjun Park
· 3_mat_mul_mgp_2026 (MGP) - Yongjun Park
댓글
댓글 쓰기