3.2 OpenMPによる並列計算

3.2.1 OpenMPによる並列計算

リスト3-2-1にSOR法のCPU版の反復計算部分を示します。
15行目でOpenMPを用いてループiを並列化しています。 またreduction句を用いて残差の和を求めています[6]。
19行目でred-black法の判定を行っています。


リスト3-2-1 OpenMPを用いた反復計算部の並列化(sor.cの一部, CPU版, vectorモード)
     1	#include "ost.h"
     2	#include "ost_prototype.h"
     3	
     4	static real_t update_rb_vector(int oe)
     5	{
     6		const real_t omega = (real_t)Solver.omega;
     7	
     8		const int64_t ni = D(1, 0, 0) - D(0, 0, 0);
     9		const int64_t nj = D(0, 1, 0) - D(0, 0, 0);
    10		const int64_t nk = D(0, 0, 1) - D(0, 0, 0);
    11	
    12		real_t res2 = 0;
    13		int64_t i;
    14	#ifdef _OPENMP
    15	#pragma omp parallel for reduction (+: res2)
    16	#endif
    17		for (        i = 0; i <= Nx; i++) {
    18		for (int64_t j = 0; j <= Ny; j++) {
    19			const int64_t kmin = ((i + j) % 2 == oe) ? 0 : 1;
    20			const int64_t kmax = Nz;
    21			const int64_t n = D(i, j, kmin);
    22			real_t *ptr = &V[n];
    23			real_t *ptr_vxp = &V[n + ni];
    24			real_t *ptr_vxm = &V[n - ni];
    25			real_t *ptr_vyp = &V[n + nj];
    26			real_t *ptr_vym = &V[n - nj];
    27			real_t *ptr_vzp = &V[n + nk];
    28			real_t *ptr_vzm = &V[n - nk];
    29			real_t *ptr_e   = &Epsr_v[n];
    30			real_t *ptr_exp = &Epsr_v[n + ni];
    31			real_t *ptr_exm = &Epsr_v[n - ni];
    32			real_t *ptr_eyp = &Epsr_v[n + nj];
    33			real_t *ptr_eym = &Epsr_v[n - nj];
    34			real_t *ptr_ezp = &Epsr_v[n + nk];
    35			real_t *ptr_ezm = &Epsr_v[n - nk];
    36			real_t rxp = RXp[i];
    37			real_t rxm = RXm[i];
    38			real_t ryp = RYp[j];
    39			real_t rym = RYm[j];
    40			real_t *ptr_rzp = &RZp[kmin];
    41			real_t *ptr_rzm = &RZm[kmin];
    42			int    *ptr_id = &idVolt_v[n];
    43			for (int64_t k = kmin; k <= kmax; k += 2) {
    44				const real_t e = *ptr_e;
    45	
    46				const real_t axp = (e + (*ptr_exp)) * rxp;
    47				const real_t axm = (e + (*ptr_exm)) * rxm;
    48				const real_t ayp = (e + (*ptr_eyp)) * ryp;
    49				const real_t aym = (e + (*ptr_eym)) * rym;
    50				const real_t azp = (e + (*ptr_ezp)) * (*ptr_rzp);
    51				const real_t azm = (e + (*ptr_ezm)) * (*ptr_rzm);
    52	
    53				const real_t res = (*ptr_id ? 0 : 1) * omega * ((
    54					axp * (*ptr_vxp) +
    55					axm * (*ptr_vxm) +
    56					ayp * (*ptr_vyp) +
    57					aym * (*ptr_vym) +
    58					azp * (*ptr_vzp) +
    59					azm * (*ptr_vzm)) / (axp + axm + ayp + aym + azp + azm) - (*ptr));
    60	
    61				*ptr += res;
    62				res2 += res * res;
    63	
    64				ptr += 2;
    65				ptr_vxp += 2;
    66				ptr_vxm += 2;
    67				ptr_vyp += 2;
    68				ptr_vym += 2;
    69				ptr_vzp += 2;
    70				ptr_vzm += 2;
    71				ptr_e += 2;
    72				ptr_exp += 2;
    73				ptr_exm += 2;
    74				ptr_eyp += 2;
    75				ptr_eym += 2;
    76				ptr_ezp += 2;
    77				ptr_ezm += 2;
    78				ptr_id += 2;
    79				ptr_rzp += 2;
    80				ptr_rzm += 2;
    81			}
    82		}
    83		}
    84	
    85		return res2;
    86	}


3.2.2 OpenMPの計算時間

表3-2-1にnovectorモードとvectorモードでスレッド数を変えたときの計算時間を示します。
これからスレッド数がある程度以上増えても速くならないことがわかります。 特にvectorモードの速度低下が顕著です。

表3-2-1 OpenMPの計算時間(FOCUSスパコンFシステム, benchmark600)
スレッド数novectorモードvectorモード
4323.0秒299.4秒
10106.2秒113.5秒
20 77.0秒152.1秒
40 96.9秒156.3秒