5선 대각 행렬을 위한 병렬 알고리즘

우 화성, 강 인석

포항 공과 대학교 화학공학과 및 기계,산업공학부

Parallel algorithm for pentadiagonal matrix

H.S. Woo, I.S. Kang

Dept. of Chem. Eng. and Division of Mech. and Indust. Eng., POSTECH

서론

최근 전산 수치 해석 분야에서는 대용량의 계산을 빠르게 처리하는 병렬 컴퓨터가 다양하게 이용되고 있다. 병렬 컴퓨터가 기존의 일반 컴퓨터의 단점인 제한적 메모리 용량과 부족한 계산 속도의 한계를 극복해주는 대안으로 제시되고 있기 때문이다. 따라서 병렬 컴퓨터는 복잡하고 많은 시간이 소요되는 시스템이나 전산 유체 역학 문제를 해석하는데 있어 가장 중요한 기술로 인식되고 있다.[1]

시스템에 관한 미분 방정식이 주어졌을 때, 이를 차분화하면 선형 미분 방정식의 형태가 되며 차분화 하는 방식에 따라서 3선 대각 행렬이나 5선 대각 행렬과 같은 띠 행렬의 형태가 된다. 2D 또는 3D와 같은 높은 차원의 시스템을 전산 묘사하거나 격자를 조밀하게 구성하여 풀게 될 경우 대량의 선형방정식들을 풀어야 함으로 많은 계산과 메모리를 요구하게 된다. 따라서 띠 행렬을 위한 효율적인 병렬 알고리즘이 필요한 것이다.

기존의 띠 행렬을 위한 병렬 알고리즘에는 우선 Factorization, Cyclic reduction, Recursive doubling과 Divide and Conquer 계열 알고리즘으로 나누어 생각 할 수 있다. Factorization 알고리즘은 수학 팩키지인 LINPACK[2]에서 사용되고 있으며 Cyclic reduction 알고리즘은 벡터 컴퓨터에 적합한 방식으로 Dubois and Rodrigue[3]에 의해 일반적인 띠 행렬에 적용 가능한 알고리즘이 개발되었다. Recursive doubling 알고리즘은 실용성이 떨어지는 것으로 알려져 있다.[4] Divide and Conquer 알고리즘은 3선 대각 행렬을 위한 Wang's partition method[5]로부터 이를 확장하여 일반적인 띠 행렬을 위한 Meier[6]의 알고리즘이 개발되었다. 이후, 3선 대각 행렬에 관한 병렬 알고리즘은 Bondeli[7]에 의해 개선되었다.

본 연구에서는 Divide and Conquer 계열의 병렬 알고리즘 중, 효율이 우수하고 확장성이 용이한 Mattor-Williams-Hewett[8]의 3선 대각 행렬을 위한 병렬 알고리즘 개념을 이용하여 higher order upwind scheme이나 biharmonic 문제 등을 풀기 위해 꼭 필요한 5선 대각 행렬을 위한 효율적인 병렬 수퍼 컴퓨터 알고리즘을 제안하고자 한다.

기본 알고리즘

다음과 같은 N×N 선형 시스템,

[eqn: A bold X = bold R]	(1)

에서	[eqn: A= BMATRIX { {b_1 }& {c_1 }& {d_1 }# {a_2}& {b_2 }& {c_2 }& {d_2}# { }& {e_3 }& {a_3 }& {b_3}& {c_3}& {d_3 }# {}# { } & { }& { ddots }& { ddots}& {ddots}& {ddots}& {ddots}# {}# { }& {}& {}& {e_N-2}& {a_N-2}& {b_N-2}& {c_N-2}& {d_N-2}# { }& {}& {}& {}& {e_N-1}& {a_N-1}& {b_N-2}& {c_N-2}# { }& {}& {}& {}& {}& {e_N}& {a_N}& {b_N} }],	(2)

[eqn: bold X = ( bold X_1{ } bold X_2 { }cdots{ } bold X_N )^T]	(3)

[eqn: bold R = ( bold R_1{ } bold R_2{ } cdots { }bold R_N )^T]	(4)

인 경우를 생각해보자. 이때, 편의상 N=PM이고, 여기서 M은 정수라고 가정하자. 여기서, P는 프로세서의 수이다. 따라서 각 하부 시스템 p에서의 하부 행렬은 다음과 같이 정의된다.

[eqn: L_p= BMATRIX { {b_1^p }& {c_1^p }& {d_1^p }# {a_2^p }& {b_2^p }& {c_2^p }& {d_2^p}# {e_3^p}&{a_3^p}&{b_3^p}&{c_3^p}&{d_3^p}# { }# { }& {ddots }& {ddots }& {ddots}& {ddots}# {}# {}&{e_M-2^p}& {a_M-2^p}& {b_M-2^p} & {c_M-2^p}&{d_M-2}# {}&{}&{e_M-1^p}& {a_M-1^p}& {b_M-1^p} & {c_M-1^p}# {}&{}&{}&{e_M^p}&{a_M^p}&{b_M^p} }],	(p=1, 2, ..., P)	(5)

같은 식으로 [eqn: bold x]와[eqn: bold r]을 정의할 수 있다. 여기서, [eqn: a_m^p == a_M(p-1)+m],, [eqn: b_m^p == b_M(p-1)+m],[eqn: c_m^p == c_M(p-1)+m],, [eqn: b_m^p == b_M(p-1)+m],, [eqn: e_m^p == e_M(p-1)+m]에 대응된다.

각각의 하부 시스템 p에서 5개의 벡터 [eqn: bold x^R ,{ } bold x ^LH_1 , { } bold x ^LH_2 , { } bold x ^UH_1 , { } bold x ^UH_2]를 다음과 같이 정의할 때,

[eqn: L_p bold x^R = [ bold r_p ],{ }L_p bold x^LH_1 = [ 0{ }0 .... -d_M-1^p{ }0 ]^T,{ }L_p bold x^LH_2 =[0{ }0...{ }0{ }-d_M^p ]^T]	(6)

[eqn: L_p bold x^UH_1 = [-e_1^p{ }0.....0{ }0]^T ,{ } L_p bold x^UH_2 = [0{ }-e_2^p ... { }0{ }0]^T]	(7)

일반해는 다음의 선형 관계식으로 정의된다.

[eqn: bold x = bold x^R + xi^LH_1 bold x ^LH_1 + xi^LH_2 bold x ^LH_2 + xi ^UH_1 bold x ^UH_1 + xi^UH_2 bold x ^UH_2]	(8)

여기서 [eqn: xi^LH_1 ,{ } xi^LH_2 ,{ } xi ^UH_1 , { } xi ^UH_2] 는 이웃한 하부 시스템의 해와의 관계를 보정 해주는 계수로써 수식 (8)을 수식 (1)에 대입하여 다음과 같은 4P-4의 계수를 결정해주는 □reduced□ 시스템을 얻을  수 있다.

[eqn: BMATRIX { {k_1^LH_1}& {k_1^LH_2 }& {-1 } # {x_1,M^LH_1 }& {x_1,M^LH_2 }& {0 }& {-1}# {}&{}&{ddots}# {-1 }& {0 }& {x_i,1^UH_1 }& {x_i,1^UH_2}& {x_i,1^LH_1}& {x_i,1^LH_2}# { }& {-1 }& {g_i^UH_1}& {g_i^UH_2}& {g_i^LH_1}& {g_i^LH_2}# { }& { }& {k_i^UH_1}& {k_i^UH_2}& {k_i^LH_1}& {k_i^LH_2}&{-1}# { }& { }& {x_i,M^UH_1}& {x_i,M^UH_2}& {x_i,M^LH_1}& {x_i,M^LH_2}& {0}& {-1}# {}&{}&{}&{}&{ddots}# { }& { }& { }&{}& {-1}& {0}& {x_P,1^UH_1}& {x_P,1^UH_2}# {}&{}&{}&{ }&{}& {-1}& {g_P^UH_1}& {g_P^UH_2}# }][eqn: BMATRIX { {xi_1^LH_1}# {xi_1^LH_2}#{ DOTSVERT }# {xi_i^UH_1 }# {xi_i^UH_2 }# {xi_i^LH_1 }# {xi_i^LH_2 }# { DOTSVERT }# {xi_P^UH_1 }# {xi_P^UH_2 } }]=[eqn: BMATRIX { {-k_1^R }# {-x_1,M^R }#{ DOTSVERT } #{-x_i,1^R }# {-g_i^R }# {-k_i^R }# {-x_i,M^R }#{dotsvert}# {-x_P,1^R }# {-g_P^R } }]	(9)

(i=2, 3, ..., P-1)

여기서,

[eqn: k_i^R = a_1^i+1 over e_1^i+1 x_i,M^R + x_i,M-1^R],[eqn: k_i^LH_1 = a_1^i+1 over e_1^i+1 x_i,M^LH_1 +x_i,M-1^LH_1 ,{ } k_i^LH_2 = a_1^i+1 over e_1^i+1 x_i,M^LH_2 +x_i,M-1^LH_2] [eqn: k_i^UH_1 = a_1^i+1 over e_1^i+1 x_i,M^UH_1 +x_i,M-1^UH_1 ,{ } k_i^UH_2 = a_1^i+1 over e_1^i+1 x_i,M^UH_2 + x_i,M-1^UH_2] 이고

[eqn: g_j^R = c_m^j-1 over d_m^j-1 x_j,1^R + x_j,2^R], [eqn: g_j^LH_1 = c_M^j-1 over d_M^j-1 x_j,1^LH_1 +x_j,2^LH_1 ,{ } g_j^LH_2 = c_M^j-1 over d_M^j-1 x_j,1^LH_2 +x_j,2^LH_2]

[eqn: g_j^UH_1 = c_M^j-1 over d_M^j-1 x_i,1^UH_1 +x_j,2^UH_1 ,{ } g_j^UH_2 = c_M^j-1 over d_M^j-1 x_j,1^UH_2 + x_j,2^UH_2] 이며,

[eqn: x_p,m]은 p번째 하부 행렬의 m번째 요소를 나타낸다.

이를 모두 수행하면 우리가 얻고자 하는 해를 (8)로부터 얻을 수 있게된다. 이상의 알고리즘의 개요를 다시 설명하자면,

(ⅰ)각각의 프로세서들은 [eqn: bold x^R ,{ } bold x ^LH_1 , { } bold x ^LH_2 , { } bold x ^UH_1 , { } bold x ^UH_2] 값들을 수식 (6), (7)로부터 동시에 계산하고, (ⅱ)수식 (9)와 같은 □reduced□ 시스템을 구축하여 [eqn: xi^LH_1 ,{ } xi^LH_2 ,] [eqn: xi^UH_1 ,{ } xi^UH_2]값들을 구한다. (ⅲ)계산된 결과들로부터 각각의 프로세서들은 (8)의 관계식으로부터 해를 구할 수 있다.

하부 행렬의 해 및 □reduced□ 시스템의 계산

수식 (6), (7)을 계산하여 [eqn: bold x^R ,{ } bold x ^LH_1 , { } bold x ^LH_2 , { } bold x ^UH_1 , { } bold x ^UH_2] 값들을 구하는데 있어, 각각의 하부 시스템에 주어진 하부 행렬은 동일하고 단지 좌변의 벡터 값들만이 다르므로 다음과 같은 효율적인 방법으로 구할 수 있다.

(1)First forward elimination

[eqn: r_p,2^UH_2 = - e_2^p],

[eqn: a_i^p = a_i^p - e_i^p b_i-1^p / a_i-1^p], [eqn: b_i^p = b_i^p - e_i^p c_i-1^p / a_i-1^p],[eqn: c_i^p = c_i^p - e_i^p d_i-1^p / a_i-1^p],

[eqn: r_i^p = r_i^p - e_i^p r_i-1^p /a_i-1^p], [eqn: r_p,i^UH_2 = -e_i^p r_p,i-1^UH_2 /a_i-1^p] (i=3,4, .... , M),

[eqn: r_p,M-1^LH_1 = - d_M-1^p], [eqn: r_p,M^LH_1 = -e_M^p d_M-1^p / a_M-1^p],

[eqn: r_p,M^LH_2 = -d_M^p].

(2)Second forward elimination

[eqn: c_1^p = c_1^p / b_1^p], [eqn: d_1^p =d_1^p /b_1^p], [eqn: r_1^p = r_1^p / b_1^p],

[eqn: r_p,1^UH_1 = -e_1^p / b_1^p], [eqn: r_p,2^UH_1 = -a_2^p r_p,1^UH_1 / (b_2^p - a_2^p c_1^p )],

[eqn: r_p,1^UH_2 = 0], [eqn: r_p,2^UH_2 = -e_2^p / (b_2^p -a_2^p c_1^p )],

[eqn: c_i^p = (c_i^p - a_i^p d_i-1^p )/ ( b_i^p -a_i^p c_i-1^p )], [eqn: d_i^p = d_i^p /(b_i^p - a_i^p - c_i-1^p )],

[eqn: r_i^p = (r_i^p - a_i^p r_i-1^p )/ (b_i^p -a_i^p c_i-1^p )] (i=2,3, ..., M),

[eqn: r_p,i^UH_1 = -a_i^p r_p,i-1^UH_1 / (b_i^p -a_i^p c_i-1^p )], [eqn: r_p,i^UH_2 = (r_p,i^UH_2 -a_i^p r_p,i-1^UH_2 ) / (b_i^p -a_i^p c_i-1^p )] (i=3,4, .. , M),

[eqn: r_p,M-1^LH_1 = r_p,M-1^LH_1 / (b_M-1^p -a_M-1^p c_M-2^p )], [eqn: r_p,M^LH_1 =(r_p,M^LH_1 -a_M^p r_p,M-1^LH_1 ) / (b_M^p - a_M^p c_M-1^p )],

[eqn: r_p,M^LH_2 = r_p,M^LH_2 / (b_M^p - a_M^p c_M-1^p )].

(3)Backward substitution

[eqn: x_p,M^R = r_M^p], [eqn: x_p,M-1^R = r_M-1^p - c_M-1^p x_p,M^R],

[eqn: x_p,M^LH_1 = r_p,M^LH_1], [eqn: x_p,M-1^LH_1 = r_p,M-1^LH_1 - c_i^p x_p,M^LH_1],

[eqn: x_p,M^LH_2 = r_p,M^LH_2], [eqn: x_p,M-1^LH_2 = -c_M-1^p x_p,M^LH_2],

[eqn: x_p,M^UH_1 = r_p,M^UH_1], [eqn: x_p,M-1^UH_1 = r_p,M-1^UH_1 - c_M-1^p x_p,M^UH_1],

[eqn: x_p,M^UH_2 = r_p,M^UH_2], [eqn: x_p,M-1^UH_2 = r_p,M-1^UH_2 - c_M-1^p x_p,M^UH_2],

[eqn: x_p,i^R = r_i^p - c_i^p x_p,i+1^R - d_i^p x_p,i+2^R],

[eqn: x_p,i^LH_1 = -c_i^p x_p,i+1^LH_1 -d_i^p x_p,i+2^LH_1], [eqn: x_p,i^LH_2 = -c_i^p x_p,i+1^LH_2 -d_i^p x_p,i+2^LH_2],

[eqn: x_p,i^UH_1 = r_p,i^UH_1 - c_i^p x_p,i+1^UH_1 - d_i^p x_p,i+2^UH_1], [eqn: x_p,i^UH_2 = r_p,i^UH_2 - c_i^p x_p,i+1^UH_2 - d_i^p x_p,i+2^UH_2] (i=M-2,M-1, ...., 1).

수식 (9)는 5선 대각 행렬 형식에서 조금 벗어나 있으나, general thomas algorithm[9]을 변형하여 쉽게 풀 수 있다.

알고리즘의 효율

5선 대각 행렬을 위한 병렬 알고리즘의 효율에 대해서 알아보기 위해서 각 계산 과정의 floating point는 다음과 같다. 여기서 M은 하부 행렬의 크기이며 P는 프로세서의 수이다.

(ⅰ) 5개의 벡터인 [eqn: bold x^R ,{ } bold x ^LH_1 , { } bold x ^LH_2 , { } bold x ^UH_1 , { } bold x ^UH_2]를 구하기 위해서, 46M

(ⅱ) "reduced" 시스템을 LU decomposition으로 계산하여 구하는 경우, 122P

(ⅲ) 최종 해를 구하기 위한 수식 (8)를 계산하기 위해서, 8M.

따라서, 하나의 floating point를 계산하기 위해 필요한 시간을 [eqn: t_b]라고 할 때, (54M+122p)[eqn: t_b]가 계산에 소요되는 시간이 된다.

그 외, "reduced" 시스템을 구성하고 구해진 계수들을 분배하기 위해서 필요한 프로세서간의 통신에 30P[eqn: t_p]+[eqn: (log_2 P) t_c]의 시간이 필요하다. 여기서 [eqn: t_c]는 하나의 통신 채널이 열리기 위해 필요한 시간을 나타내며, [eqn: t_p]는 열린 통신 채널을 통해서 단일 정보를 보내기 위해 필요한 시간을 나타낸다. 따라서, 일반적인 LU decomposition 알고리즘에 대해서 본 연구의 병렬 알고리즘의 효율은 다음과 같다.

[eqn: epsilon _p = 25 over { 54 + 122P^2 / N P log_2 P t_c / Nt_b + 30P^2 t_p / Nt_b }]	(10)

결론

기존의 3선 대각 행렬을 위한 다양한 알고리즘 중 Divide and Conquer 계열인 Mattor-Williams-Hewett의 병렬 알고리즘[8]의 개념을 확장하여 본 연구에서는 5선 대각 행렬을 위한 효율적인 병렬 알고리즘을 제안하였다. (10)에서 얻은 효율 결과에 따르면, 본 연구의 알고리즘에 관한 효율은 Mattor- Williams-Hewett의 3선 대각 행렬을 위한 병렬 알고리즘[8]의 효율과 거의 같다. 또한, 주어진 행렬의 크기가 매우 크고 프로세서의 수가 충분히 작거나 서로간의 통신속도가 매우 빠른 경우 [eqn: epsilon _p = 25 / 54 = 0.46]이 된다. 따라서, 본 연구의 알고리즘을 이용하여 higher order upwind scheme이나 biharmonic 문제와 같이 5선 대각 행렬을 반드시 풀어야 하며 많은 양의 계산이 필요한 분야에서 유용하게 이용할 수 있겠다.

감사

아낌없는 지도와 많은 조언을 해주신 포항공대 물리학과 김재삼 교수님과 포항공대-휴렛팩커드 수퍼 컴퓨팅 센터에 감사드립니다.

참고 문헌

1. S. Ragsdale, "Parallel Programming", McGraw-Hill (1991)

2. J.J. Dongarra, J.R. Bunch, C.B. Moler and G.W. Stewart, "LINPACK User's Guide", SIAM (1979)

3. P. Dubois and G. Rodrigue, in: "High Speed Computer and Algorithm Oragnization", Academic Press (1977)

4. I.S. Duff and H.A. van der Vorst, Parallel Computing, 25, 1931 (1999)

5. H.H. Wang, ACM Trans. on Math. soft., 7, 170 (1981)

6. U. Meier, Parallel Computing, 2,33 (1985)

7. S. Bondeli, Parallel Computing, 17, 419 (1991)

8. N. Mattor, T.J. Williams and D.W. Hewett, Parallel Computing, 21, 1769 (1995)

9. C.A.J. Fletcher, "Computational Techniques for Fluid Dynamics", vol. 1, Springer-Verlag (1988)