병렬 코드를 쓰려면 일단 serial 코드를 써보는 게 좋다. 그런데 잠깐 여담으로, 병렬의 counterpart인 serial을 얘기할 때, 나는 주로 순차 또는 직렬이라고 안하고 serial이라고 하게 된다. 사실 serial computing의 의미를 생각하면 순차계산이라고 하는 것이 맞다. 그런데 이거는 병렬이라는 단어랑 라임이 별로 안맞아서 counterpart로 쓰기에는 별로 마음에 안든다. 그리고 간혹가다가 직렬계산이라고 하는 사람들도 있는데, serial computing이 하나 계산하고 그 다음에 하나 계산하고 이런 의미를 갖고 있어서 단순히 일렬로 연결되어있다는 의미의 직렬을 쓰기에는 의미상 거부감이 든다. 그렇다고 시리얼/씨리얼 이라고 한국말로 쓰면 좀 콘푸라이트 생각도 나고 단어의 생김새가 좀 없어보여서 거부감이 든다. 그래서 나는 그냥 serial이라고 하는 것 같다. 수업 들을 때도 필기할 때 노트에다가 serial이라고 영어로 적게 되는 경우가 많다.
아무튼, Jacobi iteration의 내용은 뭐 워낙 쉬우니까 생략하고 바로 serial 프로그램을 써보았다.
#include <iostream>
#include <vector>
#include <cmath>
#include <random>
#include <ctime>
#include <string>
using std::vector;
using std::cout;
int main(int argc, char* argv[]){
int N = 3;
if (argc > 1){
N = std::stoi(argv[1]);
}
vector<vector<double>> A(N, vector<double>(N));
vector<double> x_new(N);
vector<double> x_old(N, 1000.0);
vector<double> b(N);
std::mt19937 gen(1);
std::uniform_real_distribution<double> dist(-0.5, 0.5);
clock_t time_make_start = std::clock();
for (int i=0; i<N; i++){
double row_sum = 0.0;
for (int j=0; j<N; j++){
if (j!=i){
A[i][j] = dist(gen);
row_sum += std::abs(A[i][j]);
}
}
A[i][i] = row_sum + 0.0001;
b[i] = dist(gen);
}
clock_t time_make_end = std::clock();
cout << "\nElapsed time for setting A and b : " <<
(double) (time_make_end - time_make_start) / CLOCKS_PER_SEC << '\n';
int iter = 0;
int max_iter = 100000;
double error = 0;
double tol = 1e-10;
clock_t time_solve_start = std::clock();
while (iter < max_iter){
for (int i=0; i<N; i++){
double sum = 0.0;
for (int j=0; j<N; j++){
if (i != j ){
sum += A[i][j] * x_old[j];
}
}
x_new[i] = (b[i] - sum)/A[i][i];
}
double error_sum = 0.0;
for (int i=0; i<N; i++){
double diff = x_new[i] - x_old[i];
error_sum += diff * diff;
}
error = std::sqrt(error_sum);
if (error < tol) break;
iter++;
x_old = x_new;
}
clock_t time_solve_end = std::clock();
cout << "Elapsed time for Jacobi iteration : " <<
(double) (time_solve_end - time_solve_start) / CLOCKS_PER_SEC << "\n\n";
cout << "No. of iterations : " << iter << '\n';
cout << "Error : " << error << '\n';
if (N<5){
cout << '\n';
for (int i = 0; i<N; i++){
if ( i == N/2)
cout << "A = ";
else
cout << " ";
cout << "[ ";
for (int j = 0; j<N; j++){
cout << A[i][j] << " ";
}
cout << "]";
if (i == N/2)
cout << " x = ";
else
cout << " ";
cout << "[ " << x_new[i] << " ]";
if (i == N/2)
cout << " b = ";
else
cout << " ";
cout << "[ " << b[i] << " ]\n";
}
cout << '\n';
}
return 0;
}
Ax=b에서 행렬 A와 b의 값을 초기화할 때, 난수 생성기를 통해서 난수로 초기화했다. 그런데 사실 수치해석을 하다 보면, 방정식을 이산화하고 수치적분 scheme을 적용해서 행렬을 구성하면 그 결과가 tri-diagonal이나 penta-diagonal처럼 nonzero의 개수가 order of N인 sparse matrix를 자주 풀게 돼서, 이 예시처럼 nonzero로 꽉차있는 행렬은 내가 자주 푸는 행렬은 아니긴 하다. 그렇지만 일단 다 값을 줬다. 그리고 Jacobi iteration이 수렴하려면 diagonal의 절대값이 그 row에 있는 다른 요소들의 절대값의 합보다 커야 해서 off-diagonal은 모두 난수로 생성하고, diagonal은 그 행의 나머지의 절대값의 합 + 0.0001로 초기화했다. 그리고 마지막에 for문으로 쓰인 cout은 행렬의 크기가 5 by 5보다 작은 경우에 터미널에 최종 행렬과 벡터를 출력하도록 쓰인 것이다. 이는 코드가 맞는지를 판단하기 위해서이다. 그래서 먼저 크기를 3 by 3으로 테스트하면 다음과 같은 결과를 얻는다.
먼저 이 값이 맞는지를 보기 위해서, 내가 코딩한 게 아닌 다른 프로그램으로 테스트하는 게 가장 빠를 것 같아서 매트랩을 통해 테스트했다.
잘 일치한다. 그러면 이제 내가 코딩한 것에서 병렬 연산 관련 코드만 추가하면 된다.
#include <iostream>
#include <vector>
#include <cmath>
#include <random>
#include <string>
#include <omp.h>
using std::vector;
using std::cout;
int main(int argc, char* argv[]){
int N = 3;
if (argc > 1){
N = std::stoi(argv[1]);
}
vector<vector<double>> A(N, vector<double>(N));
vector<double> x_new(N);
vector<double> x_old(N, 1000.0);
vector<double> b(N);
std::mt19937 gen(1);
std::uniform_real_distribution<double> dist(-0.5, 0.5);
double time_make_start = omp_get_wtime();
#pragma omp parallel
{
#pragma omp master
cout << "Using " << omp_get_num_threads() << " threads for setting A and b\n";
#pragma omp for
for (int i=0; i<N; i++){
double row_sum = 0.0;
for (int j=0; j<N; j++){
if (j!=i){
A[i][j] = dist(gen);
row_sum += std::abs(A[i][j]);
}
}
A[i][i] = row_sum + 0.0001;
b[i] = dist(gen);
}
}
double time_make_end = omp_get_wtime();
cout << "\nElapsed time for setting A and b : " << time_make_end - time_make_start << '\n';
int iter = 0;
int max_iter = 100000;
double error = 0;
double tol = 1e-10;
double time_solve_start = omp_get_wtime();
while (iter < max_iter){
#pragma omp parallel for
for (int i=0; i<N; i++){
double sum = 0.0;
for (int j=0; j<N; j++){
if (i != j ){
sum += A[i][j] * x_old[j];
}
}
x_new[i] = (b[i] - sum)/A[i][i];
}
double error_sum = 0.0;
for (int i=0; i<N; i++){
double diff = x_new[i] - x_old[i];
error_sum += diff * diff;
}
error = std::sqrt(error_sum);
if (error < tol) break;
iter++;
x_old = x_new;
}
double time_solve_end = omp_get_wtime();
cout << "Elapsed time for Jacobi iteration : " << time_solve_end - time_solve_start << "\n\n";
cout << "No. of iterations : " << iter << '\n';
cout << "Error : " << error << '\n';
if (N<5){
cout << '\n';
for (int i = 0; i<N; i++){
if ( i == N/2)
cout << "A = ";
else
cout << " ";
cout << "[ ";
for (int j = 0; j<N; j++){
cout << A[i][j] << " ";
}
cout << "]";
if (i == N/2)
cout << " x = ";
else
cout << " ";
cout << "[ " << x_new[i] << " ]";
if (i == N/2)
cout << " b = ";
else
cout << " ";
cout << "[ " << b[i] << " ]\n";
}
cout << '\n';
}
return 0;
}이렇게 일단 코딩을 했다. 난수 배정 부분 말고 실제 Jacobi를 하는 부분에서 Jacobi iteration에 대한 반복마다 새로운 parallel을 생성해서 fork와 join으로 인한 overhead가 조금 마음에 걸리긴 하지만, 일단 했다. 그리고 serial과 parallel의 시간 차이를 비교해봤다.
Jacobi를 실제로 실행하는 부분에 대한 시간은 당연히 크게 줄었지만, A와 b에 난수를 배정하는 시간이 오히려 늘었다. 이건 생각지도 못한 부분이고 이 이유를 파악하면서 꽤나 공부가 됐다. mt19937은 난수를 생성할 때, 먼저 시드에 의해서 내부적으로 현재 상태가 결정되고, 그 상태에서 연산을 통해 새로운 값을 생성해서 상태를 갱신하고 결과를 반환한다고 한다. 나도 이 정도만 읽으면 대충만 이해하고 정확히 내부적으로 어떤 일이 일어나는지는 전혀 알 수 없다. 그렇지만, 어떤 값에 접근했다가 수정한다는 것만 이해하면 된다. 이 과정을 serial로 할 때는 문제가 안되는데, 병렬로 하면, 스레드 1번도 수정하고 스레드 2번도 수정하고 스레드 3번도 수정하고 모든 스레드가 값으 수정하려고 하는데 메모리에 있는 같은 상태를 각 캐시로 가져가서 읽었다가 쓰는 순간, cache coherence를 위해서 메모리 동기화 비용이 발생한다. 그래서 이를 해결할 수 있는 방법으로 생각한 것이, 난수 생성기를 parallel 블럭 안에서 생성해서 private variable로 사용하면 cache coherence를 유지할 필요가 없어진다.
// std::mt19937 gen(1);
// std::uniform_real_distribution<double> dist(-0.5, 0.5);
double time_make_start = omp_get_wtime();
#pragma omp parallel
{
std::mt19937 gen(1);
std::uniform_real_distribution<double> dist(-0.5, 0.5);
#pragma omp master
cout << "Using " << omp_get_num_threads() << " threads for setting A and b\n";
#pragma omp for해당 부분만 다음과 같이 바꾸면 된다. 원래 parallel 블럭 전에 있던 난수 생성기를 없애고 parallel이 시작되고 나서 난수 생성기를 각 스레드마다 생성하면 된다.
이렇게 하니까 난수로 행렬 값을 채우는 시간도 훨씬 짧아졌다. 그런데 이렇게 하니까 문제가,
생각해보니까 모든 난수 생성기가 같은 시드를 사용하기 때문에 똑같은 값이 순서대로 생성될 수 밖에 없다. 코드에서는 off-diagonal 항을 먼저 생성하고, 그 row에 있는 off-diagonal 들의 절대값의 합을 이용해서 diagonal을 정하기 때문에, off-diagonal과 b가 정확히 0.497185, 0.432557, -0.371876, 0.499041로 정확히 똑같은 값이 정확히 똑같은 순서로 생성되었고, diagonal도 정확히 같은 값이 있는 걸 볼 수 있다. 이 경우에도 Jacobi iteration이 해를 구할 수는 있지만, 좀 일반적인 경우라고 보기는 어렵다. 그래서 스레드마다 난수 생성기에 다른 시드를 줘야 하는데, 스레드마다 다른 값을 가진 private 변수를 쓰면 되고, 간단하게 스레드의 번호를 이용하면 된다고 생각했다.
#pragma omp parallel
{
int tid = omp_get_thread_num();
int seed = 100*(tid * tid + 6 * tid);
std::mt19937 gen(seed);
std::uniform_real_distribution<double> dist(-0.5, 0.5);
#pragma omp masterseed에 해당하는 표현은 그냥 내가 스레드의 아이디 가지고 아무 거나 만든 거다. parallel 시작 블럭과 master 블럭 사이의 난수 생성기 선언 부분만 위와 같이 수정했다.
난수 생성 관련 문제도 해결했고 속도도 serial보다 빨라졌다.
'C++' 카테고리의 다른 글
| 병렬계산 3. OpenMP2 (0) | 2026.03.31 |
|---|---|
| 병렬계산 2. OpenMP1 (0) | 2026.03.30 |
| 병렬계산 1. Basics (0) | 2026.03.30 |