1 クリロフ部分空間法

　本ライブラリルーチンで用いた各反復法アルゴリズムについて紹介する. 以下では,n$\times$nの実行列Aを係数行列に持つ連立一次方程式 $Ax = b$ を解くものとする.ここでは,適当な前処理行列を$K$,初期値ベクトルを$x_0$とする.

1.1 前処理付共役勾配法[1]

　共役勾配法(CG法)は対称行列向けのクリロフ部分空間法である.

• アルゴリズム

1. $\mathbf{{Compute\ }}\mathbf{r}^{\mathbf{0}}\mathbf{= b - A}\mathbf{x}^{\mathbf{0}}\mathbf{{\ for\ some\ initial\ guess\ }}\mathbf{x}^{\mathbf{- 1}}\mathbf{=}\mathbf{x}^{\mathbf{0}}\mathbf{\ }$
2. $\mathbf{p}^{\mathbf{- 1}}\mathbf{= 0}$
3. $\mathbf{\alpha}_{\mathbf{- 1}}\mathbf{= 0\ }\mathbf{\ }$
4. $\mathbf{\beta}_{\mathbf{- 1}}\mathbf{= 0}$
5. $\mathbf{{solve\ s\ from\ }}\mathbf{K}\mathbf{s = \ }\mathbf{r}^{\mathbf{0}}$
6. $\mathbf{\rho}_{\mathbf{0}}\mathbf{=}\left\langle \mathbf{s,}\mathbf{r}^{\mathbf{0}} \right\rangle$
7. $\mathbf{for\ i = 0,1,2,3}\mathbf{\ldots}$
8. 　　$\mathbf{p}^{\mathbf{i}}\mathbf{= s +}\mathbf{\beta}_{\mathbf{i - 1}}\mathbf{p}^{\mathbf{i - 1}}$
9. 　　$\mathbf{q}^{\mathbf{i}}\mathbf{= A}\mathbf{p}^{\mathbf{i}}$
10. 　　$\mathbf{\gamma =}\left\langle \mathbf{p}^{\mathbf{i}}\mathbf{,}\mathbf{q}^{\mathbf{i}} \right\rangle$
11. 　　$\mathbf{x}^{\mathbf{i}}\mathbf{=}\mathbf{x}^{\mathbf{i - 1}}\mathbf{+}\mathbf{\alpha}_{\mathbf{i - 1}}\mathbf{p}^{\mathbf{i - 1}}$
12. 　　$\mathbf{\alpha}_{\mathbf{i}}\mathbf{=}\mathbf{\rho}_{\mathbf{i}}\mathbf{\ /\ \gamma}$
13. 　　$\mathbf{r}^{\mathbf{i + 1}}\mathbf{=}\mathbf{r}^{\mathbf{i}}\mathbf{-}\mathbf{\alpha}_{\mathbf{i}}\mathbf{q}^{\mathbf{i}}$
14. 　　$\mathbf{{solve\ s\ from\ }}\mathbf{K}\mathbf{s = \ }\mathbf{r}^{\mathbf{i + 1}}$
15. 　　$\mathbf{\rho}_{\mathbf{i + 1}}\mathbf{=}\left\langle \mathbf{s,}\mathbf{r}^{\mathbf{i + 1}} \right\rangle$
16. 　　$\mathbf{{if\ }}\left\| \mathbf{r}^{\mathbf{i + 1}} \right\|\mathbf{\ /\ }\left\| \mathbf{b} \right\|\mathbf{{\ small\ enough\ the}}\mathbf{n}$
17. 　　　　$\mathbf{x}^{\mathbf{i}}\mathbf{=}\mathbf{x}^{\mathbf{i - 1}}\mathbf{+}\mathbf{\alpha}_{\mathbf{i - 1}}\mathbf{p}^{\mathbf{i - 1}}$
18. 　　　　$\mathbf{{qui}}\mathbf{t}$
19. 　　$\mathbf{{endi}}\mathbf{f}$
20. 　　$\mathbf{\beta}_{\mathbf{i}}\mathbf{=}\mathbf{\rho}_{\mathbf{i + 1}}\mathbf{\ /\ }\mathbf{\rho}_{\mathbf{i}}$
21. $\mathbf{{\ en}}\mathbf{d}$

1.2 前処理付安定化双共役勾配法[1]

　安定化双共役勾配法(Bi-CGSTAB)は非対称行列向けのクリロフ部分空間法である.

• アルゴリズム

1. $\mathbf{{Compute\ }}\mathbf{r}_{\mathbf{0}}\mathbf{= b - A}\mathbf{x}_{\mathbf{0}}\mathbf{{\ for\ some\ initial\ guess\ }}\mathbf{x}_{\mathbf{- 1}}\mathbf{=}\mathbf{x}_{\mathbf{0}}$
2. $\mathbf{p}_{\mathbf{0}}\mathbf{=}\mathbf{r}_{\mathbf{0}}$
3. $\mathbf{c}_{\mathbf{1}}\mathbf{=}\left\langle \mathbf{r}_{\mathbf{0}}\mathbf{,}\mathbf{r}_{\mathbf{0}} \right\rangle$
4. $\mathbf{for\ i = 0,1,2,3}\mathbf{\ldots}$
5. 　　$\mathbf{{solve\ }}\widehat{\mathbf{p}}\mathbf{{\ from\ K}}\widehat{\mathbf{p}}\mathbf{= \ }\mathbf{p}_{\mathbf{i}}$
6. 　　$\mathbf{q = A}\widehat{\mathbf{p}}$
7. 　　$\mathbf{c}_{\mathbf{2}}\mathbf{=}\left\langle \mathbf{r}_{\mathbf{0}}\mathbf{,q} \right\rangle$
8. 　　$\mathbf{\alpha =}\mathbf{c}_{\mathbf{1}}\mathbf{\ /\ }\mathbf{c}_{\mathbf{2}}$
9. 　　$\mathbf{e =}\mathbf{r}_{\mathbf{i}}\mathbf{- \alpha q}$
10. 　　$\mathbf{{solve\ }}\widehat{\mathbf{e}}\mathbf{{\ from\ K}}\widehat{\mathbf{e}}\mathbf{= e}$
11. 　　$\mathbf{v = A}\widehat{\mathbf{e}}$
12. 　　$\mathbf{c}_{\mathbf{3}}\mathbf{=}\left\langle \mathbf{e,v} \right\rangle\mathbf{\ /\ }\left\langle \mathbf{v,v} \right\rangle$
13. 　　$\mathbf{x}_{\mathbf{i + 1}}\mathbf{=}\mathbf{x}_{\mathbf{i}}\mathbf{+ \alpha}\widehat{\mathbf{p}}\mathbf{+}\mathbf{c}_{\mathbf{3}}\widehat{\mathbf{e}}$
14. 　　$\mathbf{r}_{\mathbf{i + 1}}\mathbf{= e -}\mathbf{c}_{\mathbf{3}}\mathbf{v}$
15. 　　$\mathbf{c}_{\mathbf{1}}\mathbf{=}\left\langle \mathbf{r}_{\mathbf{0}}\mathbf{,}\mathbf{r}_{\mathbf{i + 1}} \right\rangle$
16. 　　$\mathbf{{if\ }}\left\| \mathbf{r}_{\mathbf{i + 1}} \right\|\mathbf{\ /\ }\left\| \mathbf{b} \right\|\mathbf{{\ small\ enough\ }}\mathbf{{the}}\mathbf{n}$
17. 　　　　$\mathbf{{quit}}$
18. 　　$\mathbf{{endif}}$
19. 　　$\mathbf{\beta =}\mathbf{c}_{\mathbf{1}}\mathbf{\ /\ }\left( \mathbf{c}_{\mathbf{2}}\mathbf{c}_{\mathbf{3}} \right)$
20. 　　$\mathbf{p}_{\mathbf{i + 1}}\mathbf{=}\mathbf{r}_{\mathbf{i + 1}}\mathbf{+ \beta}\left( \mathbf{p}_{\mathbf{i}}\mathbf{-}\mathbf{c}_{\mathbf{3}}\mathbf{q} \right)$
21. $\mathbf{{end}}$

1.3 前処理付一般化最小残差法[1]

　一般化最小残差法(GMRES(m))は非対称行列向けのクリロフ部分空間法である.リスタート周期m回の反復計算で解が収束しない場合はm回の反復で得られた解を初期値として反復を再開する.PARCELでは直交化手法として古典グラムシュミット法と修正グラムシュミット法が整備されている.

• アルゴリズム

$\mathbf{H}_{\mathbf{n}}$は$\mathbf{h}_{\mathbf{j,k}}$を要素に持つ上三角行列,$\mathbf{\ }\mathbf{e}_{\mathbf{i}}$はベクトル$\mathbf{e}$の最初のi成分からなるベクトルである.

1. $\mathbf{for\ j = 0,1,2,3}\mathbf{\ldots}$
2. 　　$\mathbf{r = b - A}\mathbf{x}_{\mathbf{0}}$
3. 　　$\mathbf{v}_{\mathbf{1}}\mathbf{= - r\ /\ }\left\| \mathbf{r} \right\|$
4. 　　$\mathbf{e =}\left( \mathbf{-}\left\| \mathbf{r} \right\|\mathbf{,0,\ldots,0} \right)^{\mathbf{T}}$
5. 　　$\mathbf{n = m}$
6. 　　$\mathbf{for\ i = 1,2,3}\mathbf{\ldots}\mathbf{m}$
7. 　 　　　$\mathbf{{solve\ }}{\widehat{\mathbf{v}}}_{\mathbf{i}}\mathbf{{\ from\ K}}{\widehat{\mathbf{v}}}_{\mathbf{i}}\mathbf{=}\mathbf{v}_{\mathbf{i}}$
8. 　　 　　$\mathbf{\omega = A}{\widehat{\mathbf{v}}}_{\mathbf{i}}$
9. 　　　 　$\mathbf{for\ k = 1,2,3}\mathbf{\ldots}\mathbf{i}$
10. 　　　　　　$\mathbf{h}_{\mathbf{k,i}}\mathbf{=}\left\langle \mathbf{\omega,}\mathbf{v}_{\mathbf{k}} \right\rangle$
11. 　　　　　　$\mathbf{\omega = \omega -}\mathbf{h}_{\mathbf{k,i}}\mathbf{v}_{\mathbf{k}}$
12. 　　　　$\mathbf{{end}}$
13. 　　　　$\mathbf{h}_{\mathbf{i + 1,i}}\mathbf{=}\left\| \mathbf{\omega} \right\|$
14. 　　　　$\mathbf{v}_{\mathbf{i + 1}}\mathbf{= \omega/}\left\| \mathbf{\omega} \right\|$
15. 　　　　$\mathbf{for\ k = 1,2,3}\mathbf{\ldots}\mathbf{i - 1}$
16. 　　　　　　$\begin{pmatrix} \mathbf{h}_{\mathbf{k,i}} \\ \mathbf{h}_{\mathbf{k + 1,i}} \\ \end{pmatrix}\mathbf{=}\begin{pmatrix} \mathbf{c}_{\mathbf{i}} & \mathbf{- s}_{\mathbf{i}} \\ \mathbf{s}_{\mathbf{i}} & \mathbf{c}_{\mathbf{i}} \\ \end{pmatrix}\begin{pmatrix} \mathbf{h}_{\mathbf{k,i}} \\ \mathbf{h}_{\mathbf{k + 1,i}} \\ \end{pmatrix}$
17. 　　　　$\mathbf{{end}}$
18. 　　　　$\mathbf{c}_{\mathbf{i}}\mathbf{=}\frac{\mathbf{1}}{\sqrt{\mathbf{1 +}\left( \frac{\mathbf{h}_{\mathbf{i + 1,i}}}{\mathbf{h}_{\mathbf{i,i}}} \right)^{\mathbf{2}}}}$
19. 　　　　$\mathbf{s}_{\mathbf{i}}\mathbf{= -}\frac{\mathbf{h}_{\mathbf{i + 1,i}}}{\mathbf{h}_{\mathbf{i,i}}}\frac{\mathbf{1}}{\sqrt{\mathbf{1 +}\left( \frac{\mathbf{h}_{\mathbf{i + 1,i}}}{\mathbf{h}_{\mathbf{i,i}}} \right)^{\mathbf{2}}}}$
20. 　　　　$\begin{pmatrix} \mathbf{e}_{\mathbf{i}} \\ \mathbf{e}_{\mathbf{i + 1}} \\ \end{pmatrix}\mathbf{=}\begin{pmatrix} \mathbf{c}_{\mathbf{i}} & \mathbf{- s}_{\mathbf{i}} \\ \mathbf{s}_{\mathbf{i}} & \mathbf{c}_{\mathbf{i}} \\ \end{pmatrix}\begin{pmatrix} \mathbf{e}_{\mathbf{i}} \\ \mathbf{e}_{\mathbf{i + 1}} \\ \end{pmatrix}$
21. 　　　　$\mathbf{h}_{\mathbf{i,i}}\mathbf{=}\mathbf{c}_{\mathbf{i}}\mathbf{h}_{\mathbf{i,i}}\mathbf{- \ }\mathbf{s}_{\mathbf{i}}\mathbf{h}_{\mathbf{i + 1,i}}$
22. 　　　　$\mathbf{h}_{\mathbf{i + 1,i}}\mathbf{= 0}$
23. 　　　　$\mathbf{{if\ }}\left\| \mathbf{e}_{\mathbf{i + 1}} \right\|\mathbf{\ /\ }\left\| \mathbf{b} \right\|\mathbf{{\ small\ enough\ }}\mathbf{{then}}$
24. 　　　　　　$\mathbf{n = i}$
25. 　　　　　　$\mathbf{{exit}}$
26. 　　　　$\mathbf{{endif}}$
27. 　　$\mathbf{{end}}$
28. 　　$\mathbf{y}_{\mathbf{n}}\mathbf{=}\mathbf{H}_{\mathbf{n}}^{\mathbf{- 1}}\mathbf{e}_{\mathbf{n}}$
29. 　　$\mathbf{{solve\ }}\widehat{\mathbf{x}}\mathbf{{\ from\ K}}\widehat{\mathbf{x}}\mathbf{=}\sum_{\mathbf{k = 1}}^{\mathbf{n}}{\mathbf{y}_{\mathbf{k}}\mathbf{v}_{\mathbf{k}}}$
30. 　　$\mathbf{x}_{\mathbf{n}}\mathbf{=}\mathbf{x}_{\mathbf{0}}\mathbf{+}\widehat{\mathbf{x}}$
31. 　　$\mathbf{x}_{\mathbf{0}}\mathbf{=}\mathbf{x}_{\mathbf{n}}$
32. $\mathbf{{end}}$

1.4 前処理付通信削減型一般化最小残差法 [2,3]

　通信削減型一般化最小残差法(CA-GMRES(s,t))は非対称行列向けのクリロフ部分空間法である. GMRES(m)のs回反復に相当する計算を1反復で行う. t回の反復で収束しない場合はt回の反復で得られた近似解を初期解として反復を再開する. m = t $\times$ sがGMRES(m)のリスタート周期に相当する. 数学的にはCA-GMRES(s,t)とGMRES(m)の収束性は同等であり,通信削減型のQR分解を適用することでGMRES(m)よりも集団通信の回数を削減する. ただし,一度にs本の基底ベクトルを生成するため,sが過大な場合には,丸め誤差によって基底ベクトルの直交性が崩れることでGMRESよりも収束性が悪くなる.PARCELでは基底ベクトルの直交性改善のために単基底に加えてNewton基底が整備されている.

• アルゴリズム

ここで,$\mathbf{e}$は単位ベクトル, ${\widetilde{\mathbf{\rho}}}_{\mathbf{k}}$は$\mathbf{R}_{\mathbf{k}}$のs行s列要素, $\mathbf{E}$は反復の過程で得られるヘッセンベルグ行列から求めた固有値である.

1. $\mathbf{B}_{\mathbf{k}}^{\mathbf{\sim}}\mathbf{= \lbrack}\mathbf{e}_{\mathbf{2}}\mathbf{,}\mathbf{e}_{\mathbf{3}}\mathbf{,\ldots,}\mathbf{e}_{\mathbf{s + 1}}\mathbf{\rbrack}$
2. $\mathbf{for\ j = 0,1,2,3\ldots}$
3. 　　$\mathbf{r = b - A}\mathbf{x}_{\mathbf{0}}$
4. 　　$\mathbf{v}_{\mathbf{1}}\mathbf{= r\ /\ }\left\| \mathbf{r} \right\|$
5. 　　$\mathbf{\zeta =}\left( \left\| \mathbf{r} \right\|\mathbf{,0,\ldots,0} \right)^{\mathbf{T}}$
6. 　　$\mathbf{for\ k = 0,1\ldots,t - 1}$
7. 　　　　$\mathbf{Fix\ basis\ conversion\ matrix\ \lbrack}\mathbf{B}_{\mathbf{k}}^{\mathbf{\sim}}\mathbf{,E\rbrack}$
8. 　　　　$\mathbf{Comupute\ \ basic\ vector\ \ \lbrack\ s,}{\acute{\mathbf{V}}}_{\mathbf{k}}^{\mathbf{\sim}}\mathbf{,}\mathbf{B}_{\mathbf{k}}^{\mathbf{\sim}}\mathbf{\ \rbrack}$
9. 　　　　$\mathbf{if(k.eq.0)}$
10. 　　　　　　$\mathbf{{QR\ decomposition\ }}\mathbf{V}_{\mathbf{0}}^{\mathbf{\sim}}\mathbf{=}\mathbf{Q}_{\mathbf{0}}^{\mathbf{\sim}}\mathbf{R}_{\mathbf{0}}^{\mathbf{\sim}}$
11. 　　　　　　$\mathfrak{Q}_{\mathbf{0}}^{\mathbf{\sim}}\mathbf{=}\mathbf{Q}_{\mathbf{0}}^{\mathbf{\sim}}$
12. 　　　　　　$\mathfrak{H}_{\mathbf{0}}^{\mathbf{\sim}}\mathbf{=}\mathbf{R}_{\mathbf{0}}^{\mathbf{\sim}}\mathbf{B}_{\mathbf{0}}^{\mathbf{\sim}}\mathbf{R}_{\mathbf{0}}^{\mathbf{- 1}}$
13. 　　　　　　$\mathcal{H}_{\mathbf{0}}\mathbf{=}\mathfrak{H}_{\mathbf{0}}^{\mathbf{\sim}}$
14. 　　　　　　$\mathbf{for\ o = 1,2\ldots,s}$
15. 　　　　　　$\mathbf{Givens\ rotation\ \lbrack o,}\mathcal{H}_{\mathbf{0}}\mathbf{,\zeta\rbrack}$
16. 　　　　　　$\mathbf{{if\ }}\left\| \mathbf{\zeta}_{\mathbf{o + 1}} \right\|\mathbf{\ /\ }\left\| \mathbf{b} \right\|\mathbf{{\ small\ enough\ then}}$
17. 　　　　　　　　$\mathbf{update\ \ solution\ \ vector\lbrack s,k,o,}{\acute{\mathbf{V}}}_{\mathbf{0}}^{\mathbf{\sim}}\mathbf{,}{\acute{\mathbf{Q}}}_{\mathbf{0}}^{\mathbf{\sim}}\mathbf{,}\mathbf{R}_{\mathbf{0}}\mathbf{,}\mathcal{H}_{\mathbf{0}}\mathbf{,\zeta,}\mathbf{x}_{\mathbf{0}}\mathbf{\ \rbrack}$
18. 　　　　　　　　$\mathbf{{quit}}$
19. 　　　　　　$\mathbf{{endif}}$
20. 　　　　$\mathbf{{else}}$
21. 　　　　　　${\acute{\mathfrak{R}}}_{\mathbf{k - 1,k}}^{\mathbf{\sim}}\mathbf{=}\left( \mathfrak{Q}_{\mathbf{k - 1}}^{\mathbf{\sim}} \right)^{\mathbf{H}}{\acute{\mathbf{V}}}_{\mathbf{k}}^{\mathbf{\sim}}$
22. 　　　　　　${\acute{\mathbf{V}}}_{\mathbf{k}}^{\mathbf{\sim}}\mathbf{=}{\acute{\mathbf{V}}}_{\mathbf{k}}^{\mathbf{\sim}}\mathbf{-}\mathfrak{Q}_{\mathbf{k - 1}}^{\mathbf{\sim}}{\acute{\mathfrak{R}}}_{\mathbf{k - 1,k}}^{\mathbf{\sim}}$
23. 　　　　　　$\mathbf{{QR\ decomposition\ }}{\acute{\mathbf{V}}}_{\mathbf{k}}^{\mathbf{\sim}}\mathbf{=}{\acute{\mathbf{Q}}}_{\mathbf{k}}^{\mathbf{\sim}}{\acute{\mathbf{R}}}_{\mathbf{k}}^{\mathbf{\sim}}$
24. 　　　　　　$\mathbf{R}_{\mathbf{k}}^{\mathbf{\sim}}\mathbf{=}\begin{pmatrix} \mathbf{R}_{\mathbf{k}} & \mathbf{z}_{\mathbf{k}} \\ \mathbf{0}_{\mathbf{1,k}} & \mathbf{\rho} \\ \end{pmatrix}$
25. 　　　　　　$\mathfrak{H}_{\mathbf{k - 1,k}}^{\mathbf{\sim}}\mathbf{= -}\mathfrak{H}_{\mathbf{k - 1}}\mathfrak{R}_{\mathbf{k - 1,k}}\mathbf{R}_{\mathbf{k}}^{\mathbf{- 1}}\mathbf{+}\mathfrak{R}_{\mathbf{k - 1,k}}^{\mathbf{\sim}}\mathbf{B}_{\mathbf{k}}^{\mathbf{\sim}}\mathbf{R}_{\mathbf{k}}^{\mathbf{- 1}}$
26. 　　　　　　$\mathbf{H}_{\mathbf{k}}\mathbf{=}\mathbf{R}_{\mathbf{k}}\mathbf{B}_{\mathbf{k}}\mathbf{R}_{\mathbf{k}}^{\mathbf{- 1}}\mathbf{+}{\widetilde{\mathbf{\rho}}}_{\mathbf{k}}^{\mathbf{- 1}}\mathbf{b}_{\mathbf{k}}\mathbf{z}_{\mathbf{k}}\mathbf{e}_{\mathbf{s}}^{\mathbf{T}}\mathbf{-}\mathbf{h}_{\mathbf{k - 1}}\mathbf{e}_{\mathbf{1}}\mathbf{e}_{\mathbf{s(k - 1)}}^{\mathbf{T}}\mathfrak{R}_{\mathbf{k - 1,k}}\mathbf{R}_{\mathbf{k}}^{\mathbf{- 1}}$
27. 　　　　　　$\mathbf{h}_{\mathbf{k}}\mathbf{=}{\widetilde{\mathbf{\rho}}}_{\mathbf{k}}^{\mathbf{- 1}}\mathbf{\rho}_{\mathbf{k}}\mathbf{b}_{\mathbf{k}}$
28. 　　　　　　$\mathfrak{H}_{\mathbf{k}}^{\mathbf{\sim}}\mathbf{=}\begin{pmatrix} \mathfrak{H}_{\mathbf{k - 1}} & \mathfrak{H}_{\mathbf{k - 1,k}}^{\mathbf{\sim}} \\ \begin{matrix} \mathbf{h}_{\mathbf{k - 1}}\mathbf{e}_{\mathbf{1}}\mathbf{e}_{\mathbf{s(k - 1)}}^{\mathbf{T}} \\ \mathbf{0}_{\mathbf{1,s(k - 1)}} \\ \end{matrix} & \begin{matrix} \mathbf{H}_{\mathbf{k}} \\ \mathbf{h}_{\mathbf{k}}\mathbf{e}_{\mathbf{s}}^{\mathbf{T}} \\ \end{matrix} \\ \end{pmatrix}$
29. 　　　　　　$\mathcal{H}_{\mathbf{k}}\mathbf{=}\begin{pmatrix} \mathfrak{H}_{\mathbf{k - 1,k}}^{\mathbf{\sim}} \\ \mathbf{H}_{\mathbf{k}} \\ \mathbf{h}_{\mathbf{k}}\mathbf{e}_{\mathbf{s}}^{\mathbf{T}} \\ \end{pmatrix}$
30. 　　　　　　$\mathbf{for\ o = 1 + sk,\ldots,s(k + 1)}$
31. 　　　　　　　　$\mathbf{Givens\ rotation\ \lbrack o,}\mathcal{H}_{\mathbf{k}}\mathbf{,\zeta\rbrack}$
32. 　　　　　　　　$\mathbf{{if\ }}\left\| \mathbf{\zeta}_{\mathbf{o + 1}} \right\|\mathbf{\ /\ }\left\| \mathbf{b} \right\|\mathbf{{\ small\ enough\ then}}$
33. 　　　　　　　　　　$\mathbf{update\ \ solution\ \ vector\lbrack s,k,o,}{\acute{\mathbf{V}}}_{\mathbf{k}}^{\mathbf{\sim}}\mathbf{,}\mathfrak{Q}_{\mathbf{k}}^{\mathbf{\sim}}\mathbf{,}{\acute{\mathbf{R}}}_{\mathbf{k}}^{\mathbf{\sim}}\mathbf{,}\mathcal{H}_{\mathbf{k}}\mathbf{,\zeta,}\mathbf{x}_{\mathbf{0}}\mathbf{\ \rbrack}$
34. 　　　　　　　　　　$\mathbf{{quit}}$
35. 　　　　　　　　$\mathbf{{endif}}$
36. 　　　　　　$\mathbf{{end}}$
37. 　　　　$\mathbf{{endif}}$
38. 　　$\mathbf{{end}}$
39. 　　$\mathbf{update\ \ solution\ \ vector\lbrack s,t - 1,st,}{\acute{\mathbf{V}}}_{\mathbf{t - 1}}^{\mathbf{\sim}}\mathbf{,}\mathfrak{Q}_{\mathbf{t - 1}}^{\mathbf{\sim}}\mathbf{,}{\acute{\mathbf{R}}}_{\mathbf{t - 1}}^{\mathbf{\sim}}\mathbf{,}\mathcal{H}_{\mathbf{t - 1}}\mathbf{,\zeta,}\mathbf{x}_{\mathbf{0}}\mathbf{\ \rbrack}$
40. $\mathbf{{end}}$

• Fix basis conversion matrix $\mathbf{\lbrack s,B,E\rbrack}$

1. $\mathbf{{if\ compute\ Newton\ Basis}}$
2. 　　$\mathbf{i = 0}$
3. 　　$\mathbf{{while}}\left( \mathbf{i \leq s - 1} \right)$
4. 　　　　$\mathbf{{if}}\left( \mathbf{i.eq.s - 1} \right)\mathbf{{then}}$
5. 　　　　　　$\mathbf{B}_{\mathbf{i,i}}\mathbf{= Re\lbrack}\mathbf{E}_{\mathbf{i}}\mathbf{\rbrack}$
6. 　　　　$\mathbf{{else}}$
7. 　　　　　　$\mathbf{{if}}\left( \mathbf{{Im}}\left\lbrack \mathbf{E}_{\mathbf{i}} \right\rbrack\mathbf{.eq.0} \right)\mathbf{{then}}$
8. 　　　　　　　　$\mathbf{B}_{\mathbf{i,i}}\mathbf{= Re\lbrack}\mathbf{E}_{\mathbf{i}}\mathbf{\rbrack}$
9. 　　　　　　$\mathbf{{else}}$
10. 　　　　　　　　$\mathbf{B}_{\mathbf{i,i}}\mathbf{= Re\lbrack}\mathbf{E}_{\mathbf{i}}\mathbf{\rbrack}$
11. 　　　　　　　　$\mathbf{B}_{\mathbf{i + 1,i + 1}}\mathbf{= Re\lbrack}\mathbf{E}_{\mathbf{i}}\mathbf{\rbrack}$
12. 　　　　　　　　$\mathbf{B}_{\mathbf{i,i + 1}}\mathbf{= -}\left( \mathbf{Im\lbrack}\mathbf{E}_{\mathbf{i}}\mathbf{\rbrack} \right)^{\mathbf{2}}$
13. 　　　　　　　　$\mathbf{i = i + 1}$
14. 　　　　　　$\mathbf{{endif}}$
15. 　　　　$\mathbf{{endif}}$
16. 　　　　$\mathbf{i = i + 1}$
17. 　　$\mathbf{{end\ \ }}$
18. $\mathbf{{end\ \ }}$

• Comupute basis vector $\mathbf{\lbrack\ s,v\ ,B\rbrack}$

$\mathbf{B}$は$\mathbf{b}_{\mathbf{k,i}}$を要素に持つ行列.

1. $\mathbf{for\ k = 1,2,3\ldots s}$
2. 　　$\mathbf{{solve\ }}{\widehat{\mathbf{v}}}_{\mathbf{k - 1}}\mathbf{{\ from\ K}}{\widehat{\mathbf{v}}}_{\mathbf{k - 1}}\mathbf{=}\mathbf{v}_{\mathbf{k - 1}}$
3. 　　$\mathbf{if(k.ne.1)then}$
4. 　　　　$\mathbf{\alpha =}\mathbf{b}_{\mathbf{k - 1,k - 1}}$
5. 　　　　$\mathbf{\beta =}\mathbf{b}_{\mathbf{k - 2,k - 1}}$
6. 　　　　$\mathbf{v}_{\mathbf{k}}\mathbf{= A}{\widehat{\mathbf{v}}}_{\mathbf{k - 1}}\mathbf{- \alpha}\mathbf{v}_{\mathbf{k - 1}}\mathbf{+ \beta}\mathbf{v}_{\mathbf{k - 2}}$
7. 　　$\mathbf{{else}}$
8. 　　　　$\mathbf{\alpha =}\mathbf{b}_{\mathbf{k - 1,k - 1}}$
9. 　　　　$\mathbf{v}_{\mathbf{k}}\mathbf{= A}{\widehat{\mathbf{v}}}_{\mathbf{k - 1}}\mathbf{- \alpha}\mathbf{v}_{\mathbf{k - 1}}$
10. 　　$\mathbf{{endif}}$
11. $\mathbf{{en}}\mathbf{d}$

• Givens rotation$\mathbf{\lbrack i,}\mathcal{H,}\mathbf{\zeta\rbrack}$

　$\mathcal{H}$は$\mathcal{h}_{\mathbf{k,i}}$を要素に持つ行列.

1. $\mathbf{for\ k = 1,2,3\ldots i - 1}$
2. 　　$\begin{pmatrix} \mathcal{h}_{\mathbf{k,i}} \\ \mathcal{h}_{\mathbf{k + 1,i}} \\ \end{pmatrix}\mathbf{=}\begin{pmatrix} \mathbf{c}_{\mathbf{i}} & \mathbf{- s}_{\mathbf{i}} \\ \mathbf{s}_{\mathbf{i}} & \mathbf{c}_{\mathbf{i}} \\ \end{pmatrix}\begin{pmatrix} \mathcal{h}_{\mathbf{k,i}} \\ \mathcal{h}_{\mathbf{k + 1,i}} \\ \end{pmatrix}$
3. $\mathbf{\text{end}}$
4. $\mathbf{c}_{\mathbf{i}}\mathbf{=}\frac{\mathbf{1}}{\sqrt{\mathbf{1 +}\left( \frac{\mathcal{h}_{\mathbf{i + 1,i}}}{\mathcal{h}_{\mathbf{i,i}}} \right)^{\mathbf{2}}}}$
5. $\mathbf{s}_{\mathbf{i}}\mathbf{= -}\frac{\mathcal{h}_{\mathbf{i + 1,i}}}{\mathcal{h}_{\mathbf{i,i}}}\frac{\mathbf{1}}{\sqrt{\mathbf{1 +}\left( \frac{\mathcal{h}_{\mathbf{i + 1,i}}}{\mathcal{h}_{\mathbf{i,i}}} \right)^{\mathbf{2}}}}$
6. $\begin{pmatrix} \mathbf{\zeta}_{\mathbf{i}} \\ \mathbf{\zeta}_{\mathbf{i + 1}} \\ \end{pmatrix}\mathbf{=}\begin{pmatrix} \mathbf{c}_{\mathbf{i}} & \mathbf{- s}_{\mathbf{i}} \\ \mathbf{s}_{\mathbf{i}} & \mathbf{c}_{\mathbf{i}} \\ \end{pmatrix}\begin{pmatrix} \mathbf{\zeta}_{\mathbf{i}} \\ \mathbf{\zeta}_{\mathbf{i + 1}} \\ \end{pmatrix}$
7. $\mathcal{h}_{\mathbf{i,i}}\mathbf{=}\mathbf{c}_{\mathbf{i}}\mathcal{h}_{\mathbf{i,i}}\mathbf{- \ }\mathbf{s}_{\mathbf{i}}\mathcal{h}_{\mathbf{i + 1,i}}$
8. $\mathcal{h}_{\mathbf{i + 1,i}}\mathbf{= 0}$

• Update solution vector $\mathbf{\lbrack s,k,n,V,Q,R,}\mathcal{H,}\mathbf{\zeta,}\mathbf{x}_{\mathbf{0}}\mathbf{\ \rbrack}$

1. $\mathbf{y =}\mathcal{H}^{\mathbf{- 1}}\mathbf{\zeta}$
2. $\mathbf{{solve\ }}\widehat{\mathbf{x}}\mathbf{{\ from\ K}}\widehat{\mathbf{x}}\mathbf{=}\sum_{\mathbf{k = 0}}^{\mathbf{s}\mathbf{m}\mathbf{- 1}}\mathbf{Q}_{\mathbf{k}}\mathbf{y}_{\mathbf{k}}\mathbf{+}\sum_{\mathbf{l =}\mathbf{{sm}}}^{\mathbf{n - 1}}\mathbf{V}_{\mathbf{l}\mathbf{- sm}}\sum_{\mathbf{k =}\mathbf{{sm}}}^{\mathbf{n - 1}}{\mathbf{R}_{\mathbf{l - sm}\mathbf{,}\mathbf{k}\mathbf{- sm}}^{\mathbf{- 1}}\mathbf{y}_{\mathbf{k}}}$
3. $\mathbf{x}_{\mathbf{0}}\mathbf{=}\mathbf{x}_{\mathbf{0}}\mathbf{+}\widehat{\mathbf{x}}$

2 前処理

　一般に反復解法では,$K \approx A$となるように前処理行列を選ぶと,収束に必要な反復回数が減少する. アルゴリズムに見られる 「${{solve\ }}\widehat{{p}}{{\ from\ K}}\widehat{{p}}{= \ }{p}$」型の計算は 「${A}\widehat{p}{= \ }{p}$ を$\widehat{p}$ について近似的に解く」ことに相当する. この近似の精度を上げれば反復回数が減少するが,逆にこの部分の計算のオーバーヘッドが増加する. 並列計算でない場合には,この近似計算に不完全LU 分解を用いることが多いが,並列計算の場合には,前処理計算の並列処理手法が必要となる. 本ライブラリルーチンでは並列前処理法としてポイントヤコビ法,ブロックヤコビ法,加法的シュワルツ法,ヤコビ法を採用している.

2.1 ポイントヤコビ前処理 [1]

　前処理行列Kとして係数行列Aの対角成分のみを用いる.不完全LU分解など他の前処理法と比べると効果は小さいが, 点ヤコビ法は並列に処理され,プロセッサ間通信も不要である.

2.2 フィルイン無し不完全LU分解 (ILU(0)) [1]

　正方行列$A$を下三角行列$L$と上三角行列$U$の積の形に分解することをLU分解と呼ぶ. LU分解は元の行列が疎行列であっても,分解後の行列の零要素位置に値が入る.これはフィルインと呼ばれる. フィルインが発生することでLU分解後の行列が密行列になるためメモリ使用量が増大する問題がある. PARCELではフィルインを考慮しない不完全LU分解(ILU(0))を前処理として採用することでメモリ使用量の増大を防いでいる.

2.3 修正不完全LU分解 (D-ILU) [1]

　係数行列$A$をその対角行列$D$,下三角行列$L$,上三角行列$U$を用いて $$A = L_{A} + D_{A} + U_{A}$$ と分解する.この$L_{A}$,$U_{A}$および対角行列$D$ ($D \neq D_{A}$)を用いて,前処理行列$K$を $$K = \left( D + L_{A} \right)D^{- 1}\left( D + U_{A} \right)$$ の形に決定する.ここでは,対角行列$D$は対角成分が一致するように決定する.対角行列$D$の具体的な計算方法を次に示す. $$D_{{ii}} = A_{{ii}} - \sum_{j < i}^{}{A_{{ij}}D_{jj}^{-1}A_{{ji}}}\ \ (i = 1,\ldots,n)$$

2.4 ブロックヤコビ前処理 [1]

　前処理行列K として係数行列A のブロック対角成分をとり,ブロック内で不完全LU 分解を用いる. 各ブロックの境界を係数行列の各PE への分割と一致させることにより,通信が不要となり,完全に並列に処理できる.

2.5 ヤコビ前処理 [1]

　ヤコビ前処理は次の反復法により近似的に解く前処理である.1反復内で依存性がないためGPU向けの前処理である. $${\hat{p}}_{k+1}={\hat{p}}_k+D^{-1}\left(p-A{\hat{p}}_k\right)$$

3 QR分解

　本ライブラリルーチンで利用可能なQR分解のアルゴリズムについて説明する.

3.1 古典グラムシュミット法[1]

　グラムシュミットの直交化によるQR分解.並列性は高いが直交性が悪い.

• アルゴリズム

1. $\mathbf{f}\mathbf{or\ i = 1,\ldots,n\ do}$
2. 　　$\mathbf{for\ k = 1,i - 1\ do}$
3. 　　　　$\mathbf{R}_{\mathbf{k,i}}\mathbf{= <}\mathbf{Q}_{\mathbf{k}}\mathbf{,}\mathbf{V}_{\mathbf{i}}\mathbf{>}$
4. 　　$\mathbf{{enddo\ }}$
5. 　　$\mathbf{for\ k = 1,i - 1\ do}$
6. 　　　　$\mathbf{V}_{\mathbf{i}}\mathbf{=}\mathbf{V}_{\mathbf{i}}\mathbf{-}\mathbf{R}_{\mathbf{k,i}}\mathbf{Q}_{\mathbf{k}}$
7. 　　　　$\mathbf{{enddo\ }}$
8. 　　$\mathbf{R}_{\mathbf{i,i}}\mathbf{= <}\mathbf{V}_{\mathbf{i}}\mathbf{,}\mathbf{V}_{\mathbf{i}}\mathbf{>}$
9. 　　$\mathbf{Q}_{\mathbf{i}}\mathbf{=}\frac{\mathbf{1}}{\left\| \mathbf{V}_{\mathbf{i}} \right\|}\mathbf{V}_{\mathbf{i}}$
10. $\mathbf{{enddo\ }}$

3.2 修正グラムシュミット法[1]

　誤差の低減を目的に古典グラムシュミット法のアルゴリズムを修正したQR分解. 古典グラムシュミット法よりも直交性が高いが集団通信の回数が多い.

• アルゴリズム

1. $\mathbf{f}\mathbf{or\ i = 1,\ldots,n\ do}$
2. 　　$\mathbf{for\ k = 1,i - 1\ do}$
3. 　　　　$\mathbf{R}_{\mathbf{k,i}}\mathbf{= <}\mathbf{Q}_{\mathbf{k}}\mathbf{,}\mathbf{V}_{\mathbf{i}}\mathbf{>}$
4. 　　　　$\mathbf{V}_{\mathbf{i}}\mathbf{=}\mathbf{V}_{\mathbf{i}}\mathbf{-}\mathbf{R}_{\mathbf{k,i}}\mathbf{Q}_{\mathbf{k}}$
5. 　　$\mathbf{{enddo\ }}$
6. 　　$\mathbf{R}_{\mathbf{i,i}}\mathbf{= <}\mathbf{V}_{\mathbf{i}}\mathbf{,}\mathbf{V}_{\mathbf{i}}\mathbf{>}$
7. 　　$\mathbf{Q}_{\mathbf{i}}\mathbf{=}\frac{\mathbf{1}}{\left\| \mathbf{V}_{\mathbf{i}} \right\|}\mathbf{V}_{\mathbf{i}}$
8. $\mathbf{{enddo\ }}$

3.3 コレスキーQR法[4]

　コレスキーQRは行列積とコレスキー分解で構成されるQR分解である.演算密度が高く1回の集団通信で計算可能であるが直交性が悪い.

• アルゴリズム

1. 　$\mathbf{B =}\mathbf{V}^{\mathbf{T}}\mathbf{V}$
2. 　$\mathbf{Cholesky\ decomposition(B)}$
3. 　$\mathbf{Q = V}\mathbf{R}^{\mathbf{- 1}}$

4 データ構造

　一般に疎行列反復解法のプログラムでは,メモリ節約のため,非零要素についての情報(値,行位置,列位置) のみを格納するデータ構造を採用する.ここでは,本ライブラリルーチンで用いている疎行列格納法について説明する.

4.1 Compressed Row Storage (CRS)形式

　CRS形式(圧縮行格納形式)は,非ゼロ要素を行方向に圧縮した後,列方向にデータを並べて格納する方法である. 非ゼロ要素の列位置・値ともに全体を1 次元配列に格納する.そのため.1次元配列中での各行の開始位置を記録したポインタテーブルを用いる. 並列処理の場合には列方向に行列を分割して格納する. 列方向のポインターの並びは昇順を想定している.

• 行列構造

$$A = \left( \begin{array}{ccc} a & b & c & 0\\ 0 & d & 0 & 0\\ e & 0 & f & g\\ 0 & h & i & j \end{array} \right)$$

1プロセスの場合

• ランク0プロセス

行列の非ゼロ要素　　::　$crsA=\{a,b,c,d,e,f,g,h,i,j\}$

ポインタテーブル　　::　$crsRow\_ptr=\{1,4,5,8,11\}$

非ゼロ要素の列番号　::　$crsCol=\{1,2,3,2,1,3,4,2,3,4\}$

2プロセスの場合

• ランク0プロセス

行列の非ゼロ要素　　::　 $crsA=\{a,b,c,d\}$

ポインタテーブル　　::　$crsRow\_ptr=\{1,4,5\}$

非ゼロ要素の列番号　::　$crsCol=\{1,2,3,2\}$

• ランク1プロセス

行列の非ゼロ要素　　::　$crsA=\{e,f,g,h,i,j\}$

ポインタテーブル　　::　$crsRow\_ptr=\{1,4,7\}$

非ゼロ要素の列番号　::　$crsCol=\{1,3,4,2,3,4\}$

5 インストール

　SYCL版PARCELの導入方法について説明する.
　ヘッダーファイル./Include/parcel_sycl.hppをインクルードしてSYCL対応コンパイラでコンパイルすれば利用可能である. 以下では SGI8600 (Intel Xeon Gold 6248R + NVIDIA V100) における例を示す.

• CPUプログラムのコンパイル

  module purge
  module load gnu/9.5.0
  module load openmpi/cur
  . ~/spack/share/spack/setup-env.sh
  spack load hipsycl
  export sycl_root=`spack location -i hipsycl`
  export llvm_root=`spack location -i llvm`
  export LD_LIBRARY_PATH=$sycl_root/lib/:$llvm_root/lib/:${LD_LIBRARY_PATH}
  export MPICXX_CXX=${compiler_base}
  export OMPI_CXX=${compiler_base}
  src_name=sample.cpp
  compiler_base="acpp"
  compiler_option=" -O3 --acpp-targets=omp -mlong-double-128  -march=native -mtune=native  -Wpass-failed=transform-warning"
  compiler=mpicxx
  ${compiler} ${compiler_option} -o a.out ${src_name}

• GPUプログラムのコンパイル

  module purge
  module load gnu/9.5.0
  module load cuda/12.6
  module load openmpi/cur
  . ~/spack/share/spack/setup-env.sh
  spack load hipsycl
  export sycl_root=`spack location -i hipsycl`
  export llvm_root=`spack location -i llvm`
  export cuda_root=/home/app/cuda/12.6
  export ACPP_CUDA_PATH=$cuda_root
  export LD_LIBRARY_PATH=$sycl_root/lib/:$llvm_root/lib/:$cuda_root/lib64:${LD_LIBRARY_PATH}
  export MPICXX_CXX=${compiler_base}
  export OMPI_CXX=${compiler_base}
  src_name=sample.cpp
  compiler_base="acpp"
  compiler_option=" -O3 --acpp-targets=cuda:sm_70 -march=native -mtune=native -v"
  compiler=mpicxx
  ${compiler} ${compiler_option} -o a.out ${src_name}

6 SYCL対応コンパイラのインストール方法

6.1 AdaptiveCpp

ソースコードからインストールする場合の詳細は下記のページを参照.

https://github.com/AdaptiveCpp/AdaptiveCpp

spackによるインストールも可能である. 以下では SGI8600 (Intel Xeon Gold 6248R + NVIDIA V100) における例を示す.

• SGI8600におけるCPU環境のインストール

  module purge
  module load gnu/9.5.0
  module load openmpi/cur
  . ~/spack/share/spack/setup-env.sh
  spack compiler find
  spack external find
  spack install llvm@18
  spack install hipsycl % llvm@18

富岳やWisteria Odysseyのようにログインノードと計算ノードでCPUアーキテクチャが異なる場合には計算ノード上でインストールを行う.

• SGI8600におけるGPU環境のインストール

  module purge
  module load gnu/9.5.0
  module load cuda/12.6
  module load openmpi/cur
  . ~/spack/share/spack/setup-env.sh
  spack compiler find
  spack external find
  spack external find cuda
  spack install llvm@18
  spack install hipsycl +cuda % llvm@18

7 SYCL版PARCELクラス

　SYCL版PARCELに整備されているクラスについて説明する.

7.1 PARCELデータ型

　SYCL版PARCELはデータ型として次のデータ型(parceltype)を利用する.

• parcel_fp16 : sycl::half型
• parcel_float : float型
• parcel_double : double型
• parcel_longdouble : long double型
• parcel_doubledouble : double-double型(4倍精度)
• PARCEL_INT : 整数型

7.2 parcel::comm

SYCL版PARCEL用MPIコミュニケータクラス.

• 構文

7.3 parcel::matrix

SYCL版PARCELの疎行列を操作するクラス.

• 構文

• メンバ関数

7.4 parcel::vector

SYCL版PARCELのベクトルを操作するクラス.

• 構文

• メンバ関数

7.5 parcel::krylov::pcg

前処理付き共役勾配法クラス.

• 構文

• メンバ関数

7.6 parcel::krylov::pbicgstab

前処理付安定化双共役勾配法クラス.

• 構文

• メンバ関数

7.7 parcel::krylov::pgmres

前処理付一般化最小残差法クラス.

• 構文

• メンバ関数

7.8 parcel::krylov::pcagmres

前処理付き通信削減型一般化最小残差法クラス.

• 構文

• メンバ関数

8 使用方法

　SYCL版PARCELのサンプルプログラムを示す.

8.1 C++

  #include <mpi.h>
 
  #include "./Include/parcel_sycl.hpp"
  #include "./parcel_default.hpp"
 
  void run_parcel_sycl_main(
              MPI_Comm mpi_comm_c, // In:MPIコミュニケータ
              int n, // In:CRS行列行数
              int nnz, // In:CRS行列非零要素数
              double* crsA, // In:CRS行列要素
              int* crsCol, // In:CRS形式列番号
              int* crsRow_ptr, // In:CRS行列ポインタテーブル
              double* vecx, // In:初期ベクトル,out:解ベクトル
              double* vecb  // In:右辺ベクトル
               ) {
 
    sycl::device device;
    sycl::property::queue::in_order property;
    sycl::queue q{device,property};
    q = sycl::default_selector{};
 
    parcel::comm comm(mpi_comm_c,&q);
    int myrank   = comm.mpi_rank();
    int nprocess = comm.mpi_size();
 
    if(myrank == 0)std::cout << q.get_device().get_info() << std::endl;
 
    parcel::matrix parcel_matrix(&q);
    PARCEL_INT nlocal = n;
    PARCEL_INT nstart = nlocal*myrank;
    PARCEL_INT nend   = nlocal*(myrank+1)-1;
 
    q.wait();
    parcel_matrix.set_crsMat(
                     &comm,nstart,nend,
                     n,nnz,
                     crsRow_ptr,
                     crsCol,
                     crsA
                     );
    q.wait();
 
    parcel::vector vx(n,&q);
    parcel::vector vb(n,&q);
    vx.setValue(vecx);
    vb.setValue(vecb);
 
 
    parcel::krylov::pcg<
      parcel_solver_calc_type,
      parcel_solver_calc_precon_type,
      parcel_solver_calc_type,
      parcel_solver_calc_type
      > solver;
    solver.set_tolerance( 1.0e-9 );
    solver.set_max_iteration( MAX_ITERATION );
 
    PARCEL_INT precon_algo = parcel_solver_precon_algo;
    switch (precon_algo) {
    case parcel::PRECON_NONE:
      // none
      solver.set_precon_algo(parcel::PRECON_NONE);
      break;
    case parcel::PRECON_POINTJACOBI:
      // point jacobi precon
      solver.set_precon_algo(parcel::PRECON_POINTJACOBI);
      break;
    case parcel::PRECON_BLOCKJACOBI:
      // block jacobi precon
      solver.set_precon_algo(parcel::PRECON_BLOCKJACOBI);
      solver.set_scaling_matrix(true);
      solver.set_scaling_vector(true);
      solver.set_blockjacobi(BLOCKJACOBI_BLOCK);
      break;
    case parcel::PRECON_JACOBI:
      // jacobi precon
      solver.set_precon_algo(parcel::PRECON_JACOBI);
      solver.set_scaling_matrix(true);
      solver.set_scaling_vector(true);
      solver.set_jacobi(JACOBI_ITER);
      break;
    }
    q.wait();
 
    solver.solver<
      parcel_solver_data_type,
      parcel_solver_data_type,
      parcel_solver_data_type,
      parcel_solver_data_type,
      parcel_solver_data_precon_mat_type,
      parcel_solver_data_precon_vec_type
      >(parcel_matrix,vx,vb);
    vx.getValue(vecx);  q.wait();
    q.wait();
 
    return ;
  }

8.2 Fortran

• Fortran用wrapper関数(C++)

  extern "C" {
    void run_parcel_sycl_f(
             int* mpi_comm_f,
             int n, int nnz,
             double* crsA,
             int* crsCol,
             int* crsRow_ptr,
             double* vecx,
             double* vecb
              ) {
      MPI_Comm mpi_comm_c = MPI_Comm_f2c(mpi_comm_f[0]);
      run_parcel_sycl_main(
             mpi_comm_c,
             n, nnz,
             crsA,
             crsCol,
             crsRow_ptr,
             vecx,
             vecb );
      return ;
    }
  }
 

• Fortranプログラム

make_poisson3d_crs_get_nnz,make_poisson3d_crsは行列分割に対応したCRS形式の行列を生成するサブルーチンとする.

  module parcel_sycl_wrapper
 
    interface
       subroutine run_parcel_sycl( &
            mpi_comm_f, &
            n,nnz, &
            crsA, &
            crsCol, &
            crsRow_ptr, &
            vecx, &
            vecb &
            ) &
            bind(c,name="run_parcel_sycl_f")
         implicit none
         integer mpi_comm_f
         integer(kind=4),value :: n
         integer(kind=4),value :: nnz
         real*8 crsA(*)
         integer(kind=4) crsCol(*)
         integer(kind=4) crsRow_ptr(*)
         real*8 vecx(*)
         real*8 vecb(*)
 
       end subroutine run_parcel_sycl
    end interface
  contains
 
  end module parcel_sycl_wrapper
 
  program main
    use mpi
    use parcel_sycl_wrapper
    implicit none
    integer(kind=4) nx,ny,nz
    integer(kind=4) nx0,ny0,nz0
 
    integer(kind=4),parameter :: gnx = 16
    integer(kind=4),parameter :: gny = 16
    integer(kind=4),parameter :: gnz = 16
 
    integer(kind=4) gn,n
    integer(kind=4) nnz
    integer(kind=4) istart
 
    real*8,allocatable :: crsA(:)
    integer(kind=4),allocatable :: crsCol(:)
    integer(kind=4),allocatable :: crsRow_ptr(:)
 
    real*8,allocatable :: vecx(:)
    real*8,allocatable :: vecb(:)
 
    integer ierr
    integer npes,myrank
 
    integer i
    integer mpi_comm_f
 
    mpi_comm_f = MPI_COMM_WORLD
 
    ! MPI Init
    call MPI_Init(ierr)
    call MPI_Comm_size(mpi_comm_f,npes,ierr)
    call MPI_Comm_rank(mpi_comm_f,myrank,ierr)
 
    nx     = gnx
    ny     = gny
    nz     = gnz/npes
 
    nx0    = 1
    ny0    = 1
    nz0    = 1 + nz*myrank
 
    gn     = gnx*gny*gnz
    n      = nx * ny * nz
    istart = 1+n*myrank
 
    call make_poisson3d_crs_get_nnz( &
         nx,ny,nz,   &
         nx0,ny0,nz0,&
         gnx,gny,gnz,&
         nnz, &
         istart)
 
    allocate(crsA(nnz))
    allocate(crsCol(nnz))
    allocate(crsRow_ptr(n+1))
    allocate(vecx(n))
    allocate(vecb(n))
 
    vecx = 0.0d0
    vecb = 1.0d0
    call make_poisson3d_crs( &
         nx,ny,nz,   &
         nx0,ny0,nz0,&
         gnx,gny,gnz,&
         crsA,crsCol,crsRow_ptr,n,nnz, &
         istart)
 
    call run_parcel_sycl( &
         mpi_comm_f, &
         n,nnz, &
         crsA, &
         crsCol, &
         crsRow_ptr, &
         vecx, &
         vecb &
         )
 
    deallocate(crsA)
    deallocate(crsCol)
    deallocate(crsRow_ptr)
    deallocate(vecx)
    deallocate(vecb)
 
    call MPI_Finalize(ierr);
 
  end program main

8.3 C言語

• C言語用wrapper関数(C++)

  extern "C" {
    void run_parcel_sycl_c(
             MPI_Comm mpi_comm_c,
             int n, int nnz,
             double* crsA,
             int* crsCol,
             int* crsRow_ptr,
             double* vecx,
             double* vecb
              ) {
      run_parcel_sycl_main(
             mpi_comm_c,
             n, nnz,
             crsA,
             crsCol,
             crsRow_ptr,
             vecx,
             vecb );
      return ;
    }
  }

• C言語プログラム

make_poisson3d_crs_get_nnz_,make_poisson3d_crs_は行列分割に対応したCRS形式の行列を生成する関数とする.

  #include <stdio.h>
  #include <stdlib.h>
  #include <mpi.h>
 
  extern void make_poisson3d_crs_get_nnz_(
                       int* nx,
                       int* ny,
                       int* nz,
                       int* nx0,
                       int* ny0,
                       int* nz0,
                       int* gnx,
                       int* gny,
                       int* gnz,
                       int* nnz,
                       int* istart);
  extern void make_poisson3d_crs_(
                   int* nx,
                   int* ny,
                   int* nz,
                   int* nx0,
                   int* ny0,
                   int* nz0,
                   int* gnx,
                   int* gny,
                   int* gnz,
                   double* crsA,
                   int* crsCol,
                   int* crsRow_ptr,
                   int* n,
                   int* nnz0,
                   int* istart);
  extern void run_parcel_sycl_c(
               MPI_Comm mpi_comm_c,
               int n, int nnz,
               double* crsA,
               int* crsCol,
               int* crsRow_ptr,
               double* vecx,
               double* vecb
            );
 
  int main(int argc, char* argv[]) {
 
    int gnx = 16;
    int gny = 16;
    int gnz = 16;
 
    MPI_Init(&argc, &argv);
    MPI_Comm mpi_comm_c = MPI_COMM_WORLD;
    int npes,myrank;
    MPI_Comm_size(MPI_COMM_WORLD,&npes);
    MPI_Comm_rank(MPI_COMM_WORLD,&myrank);
 
    int nx     = gnx;
    int ny     = gny;
    int nz     = gnz/npes;
 
    int nx0    = 1;
    int ny0    = 1 ;
    int nz0    = 1 + nz*myrank;
 
    int gn     = gnx*gny*gnz;
    int n      = nx * ny * nz;
    int istart = 1+n*myrank;
    int nnz;
 
    make_poisson3d_crs_get_nnz_(
         &nx,&ny,&nz,
         &nx0,&ny0,&nz0,
         &gnx,&gny,&gnz,
         &nnz,
         &istart);
 
    double* crsA = (double*)malloc(sizeof(double)*nnz);
    int* crsCol = (int*)malloc(sizeof(int)*nnz);
    int* crsRow_ptr = (int*)malloc(sizeof(int)*(n+1));
    double* vecx = (double*)malloc(sizeof(double)*n);
    double* vecb = (double*)malloc(sizeof(double)*n);
 
    make_poisson3d_crs_(
             &nx,&ny,&nz,
             &nx0,&ny0,&nz0,
             &gnx,&gny,&gnz,
             crsA,crsCol,crsRow_ptr,
             &n,&nnz,
             &istart);
    for(int i = 0;i<n;i++){
      vecx[i] = 0.0;
      vecb[i] = 1.0;
    }
 
    run_parcel_sycl_c(
            mpi_comm_c,
            n,nnz,
            crsA,
            crsCol,
            crsRow_ptr,
            vecx,
            vecb
             );
 
    free(crsA);
    free(crsCol);
    free(crsRow_ptr);
    free(vecx);
    free(vecb);
 
    MPI_Finalize();
 
    return 0;
  }

9 参考文献

[1]R. Barret, M. Berry, T. F. Chan, et al., “Templates for the Solution of Linear Systems: Building Blocks for Iterative Methods”, SIAM（1994）
[2]M. Hoemmen, “Communication-avoiding Krylov subspace methods”. Ph.D.thesis, University of California, Berkeley（2010）
[3]Y. Idomura, T. Ina, A. Mayumi, et al., “Application of a communication-avoiding generalized minimal residual method to a gyrokinetic five dimensional eulerian code on many core platforms”,ScalA 17:8th Workshop on Latest Advances in Scalable Algorithms for Large Scale Systems, pp.1-8, (2017).
[4] A. Stathopoulos, K. Wu, “A block orthogonalization procedure with constant synchronization requirements”. SIAM J. Sci. Comput. 23, 2165–2182 (2002)