1 クリロフ部分空間法

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

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

　共役勾配法(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,2]

　安定化双共役勾配法(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,2]

　一般化最小残差法(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 前処理付通信削減型一般化最小残差法 [3,4]

　通信削減型一般化最小残差法(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}}\)

1.5 前処理付チェビシェフ基底共役勾配法 [5,6]

　チェビシェフ基底共役勾配法(CBCG)は対称行列向けのクリロフ部分空間法である.CG法のk反復に相当する計算を1反復で行うことで集団通信の回数を削減する.チェビシェフ基底を用いて基底ベクトルを生成するため行列の最大・最小固有値を必要とする.PARCELでは固有値計算法としてべき乗法と通信削減型アーノルディ法が整備されている.

• アルゴリズム

1. \(\mathbf{{Compute\ }}\mathbf{r}_{\mathbf{0}}\mathbf{= b - A}\mathbf{x}_{\mathbf{0}}\mathbf{{\ for\ some\ initial\ gue}}\mathbf{{ss\ }}\mathbf{x}_{\mathbf{- 1}}\mathbf{=}\mathbf{x}_{\mathbf{0}}\)
2. \(\mathbf{Compute\ Chebyshev\ basis\ \lbrack}\mathbf{S}_{\mathbf{0}}\mathbf{,}{\mathbf{A}\mathbf{S}}_{\mathbf{0}}\mathbf{,}\mathbf{r}_{\mathbf{0}}\mathbf{,}\mathbf{\lambda}_{\mathbf{\max}}\mathbf{,}\mathbf{\lambda}_{\mathbf{\min}}\mathbf{\rbrack}\)
3. \(\mathbf{Q}_{\mathbf{0}}\mathbf{=}\mathbf{S}_{\mathbf{0}}\)
4. \(\mathbf{for\ i = 0,1,2,3\ldots}\)
5. \(\mathbf{{Compute\ }}\mathbf{Q}_{\mathbf{i}}^{\mathbf{T}}\mathbf{A}\mathbf{Q}_{\mathbf{i}}\mathbf{\ ,\ }\mathbf{Q}_{\mathbf{i}}^{\mathbf{T}}\mathbf{r}_{\mathbf{{ik}}}\)
6. 　　\(\mathbf{\alpha}_{\mathbf{i}}\mathbf{= \ }\left( \mathbf{Q}_{\mathbf{i}}^{\mathbf{T}}\mathbf{A}\mathbf{Q}_{\mathbf{i}} \right)^{\mathbf{- 1}}\left( \mathbf{Q}_{\mathbf{i}}^{\mathbf{T}}\mathbf{r}_{\mathbf{{ik}}} \right)\)
7. 　　\(\mathbf{x}_{\left( \mathbf{i + 1} \right)\mathbf{k}}\mathbf{=}\mathbf{x}_{\mathbf{{ik}}}\mathbf{+}\mathbf{Q}_{\mathbf{i}}\mathbf{\alpha}_{\mathbf{i}}\)
8. 　　\(\mathbf{r}_{\left( \mathbf{i + 1} \right)\mathbf{k}}\mathbf{=}\mathbf{r}_{\mathbf{{ik}}}\mathbf{- A}\mathbf{Q}_{\mathbf{i}}\mathbf{\alpha}_{\mathbf{i}}\)
9. 　　\(\mathbf{{if\ }}\left\| \mathbf{r}_{\left( \mathbf{i + 1} \right)\mathbf{k}} \right\|\mathbf{\ /\ }\left\| \mathbf{b} \right\|\mathbf{{\ small\ enough\ }}\mathbf{{then}}\)
10. 　　\(\mathbf{Compute\ Chebyshev\ basis\ \lbrack}\mathbf{S}_{\mathbf{i + 1}}\mathbf{,}{\mathbf{A}\mathbf{S}}_{\mathbf{i + 1}}\mathbf{,}\mathbf{r}_{\left( \mathbf{i + 1} \right)\mathbf{k}}\mathbf{,}\mathbf{\lambda}_{\mathbf{\max}}\mathbf{,}\mathbf{\lambda}_{\mathbf{\min}}\mathbf{\rbrack}\)
11. 　　\(\mathbf{{Compute\ }}\mathbf{Q}_{\mathbf{i}}^{\mathbf{T}}\mathbf{A}\mathbf{S}_{\mathbf{i + 1}}\)
12. 　　\(\mathbf{B}_{\mathbf{i}}\mathbf{=}\left( \mathbf{Q}_{\mathbf{i}}^{\mathbf{T}}\mathbf{A}\mathbf{Q}_{\mathbf{i}} \right)^{\mathbf{- 1}}\left( \mathbf{Q}_{\mathbf{i}}^{\mathbf{T}}\mathbf{A}\mathbf{S}_{\mathbf{i + 1}} \right)\mathbf{\ }\)
13. 　　\(\mathbf{Q}_{\mathbf{i + 1}}\mathbf{=}\mathbf{S}_{\mathbf{i + 1}}\mathbf{-}\mathbf{Q}_{\mathbf{i}}\mathbf{B}_{\mathbf{i}}\)
14. 　　\(\mathbf{A}\mathbf{Q}_{\mathbf{i + 1}}\mathbf{=}\mathbf{{AS}}_{\mathbf{i + 1}}\mathbf{-}\mathbf{{AQ}}_{\mathbf{i}}\mathbf{B}_{\mathbf{i}}\)
15. \(\mathbf{{end}}\)

• Compute Chebyshev basis \(\mathbf{\lbrack S,AS,r,}\mathbf{\lambda}_{\mathbf{\max}}\mathbf{,}\mathbf{\lambda}_{\mathbf{\min}}\mathbf{\rbrack}\)

\(\mathbf{\lambda}_{\mathbf{\max}}\),\(\mathbf{\lambda}_{\mathbf{\min}}\)は\(\mathbf{AK}^{-1}\)の最大及び最小固有値.

1. \(\mathbf{\eta = 2\ /\ (}\mathbf{\lambda}_{\mathbf{\max}}\mathbf{-}\mathbf{\lambda}_{\mathbf{\min}}\mathbf{)}\)
2. \(\mathbf{\zeta =}\left( \mathbf{\lambda}_{\mathbf{\max}}\mathbf{+}\mathbf{\lambda}_{\mathbf{\min}} \right)\mathbf{\ /\ (}\mathbf{\lambda}_{\mathbf{\max}}\mathbf{-}\mathbf{\lambda}_{\mathbf{\min}}\mathbf{)}\)
3. \(\mathbf{s}_{\mathbf{0}}\mathbf{=}\mathbf{r}_{\mathbf{0}}\mathbf{\ }\)
4. \(\mathbf{{solve\ }}{\widetilde{\mathbf{s}}}_{\mathbf{0}}\mathbf{{\ from\ K}}{\widetilde{\mathbf{s}}}_{\mathbf{0}}\mathbf{= \ }\mathbf{s}_{\mathbf{0}}\)
5. \(\mathbf{s}_{\mathbf{1}}\mathbf{= \eta A}{\widetilde{\mathbf{s}}}_{\mathbf{0}}\mathbf{- \zeta}\mathbf{s}_{\mathbf{0}}\)
6. \(\mathbf{{solve\ }}{\widetilde{\mathbf{s}}}_{\mathbf{1}}\mathbf{{\ from\ K}}{\widetilde{\mathbf{s}}}_{\mathbf{1}}\mathbf{= \ }\mathbf{s}_{\mathbf{1}}\)
7. \(\mathbf{for\ j = 2,3,\ldots,k}\)
8. 　　\(\mathbf{s}_{\mathbf{j}}\mathbf{= 2\eta A}{\widetilde{\mathbf{s}}}_{\mathbf{j}}\mathbf{- 2\zeta}\mathbf{s}_{\mathbf{j - 1}}\mathbf{-}\mathbf{s}_{\mathbf{j - 2}}\)
9. 　　\(\mathbf{{solve\ }}{\widetilde{\mathbf{s}}}_{\mathbf{j}}\mathbf{{\ from\ K}}{\widetilde{\mathbf{s}}}_{\mathbf{j}}\mathbf{= \ }\mathbf{s}_{\mathbf{j}}\)
10. \(\mathbf{{end}}\)
11. \(\mathbf{S = (}{\widetilde{\mathbf{s}}}_{\mathbf{0}}\mathbf{,}{\widetilde{\mathbf{s}}}_{\mathbf{1}}\mathbf{,\ldots,}{\widetilde{\mathbf{s}}}_{\mathbf{k - 1}}\mathbf{)}\)
12. \(\mathbf{AS = (A}{\widetilde{\mathbf{s}}}_{\mathbf{0}}\mathbf{,}{\mathbf{A}\widetilde{\mathbf{s}}}_{\mathbf{1}}\mathbf{,\ldots,}{\mathbf{A}\widetilde{\mathbf{s}}}_{\mathbf{k - 1}}\mathbf{)}\)

1.6 通信削減型アーノルディ法 [3]

　通信削減型アーノルディ法(CA-Arnoldi(s,t))は非対称行列向けの固有値計算アルゴリズムである. アーノルディ法のs回反復に相当する計算を1反復で行う.t回反復することでt \(\times\) sのヘッセンベルグ行列を作成して固有値と固有ベクトルを計算する. 通信削減型のQR分解を適用することでアーノルディ法よりも集団通信の回数を削減する.

• アルゴリズム

ここで,\(\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{r = }\mathbf{x}_{\mathbf{0}}\)
3. 　　\(\mathbf{\zeta =}\left( \left\| \mathbf{r} \right\|\mathbf{,0,\ldots,0} \right)^{\mathbf{T}}\)
4. 　　\(\mathbf{for\ k = 0,1\ldots,t - 1}\)
5. 　　　　\(\mathbf{Fix\ basis\ conversion\ matrix\ \lbrack}\mathbf{B}_{\mathbf{k}}^{\mathbf{\sim}}\mathbf{,E\rbrack}\)
6. 　　　　\(\mathbf{Comupute\ \ basic\ vector\ \ \lbrack\ s,}{\acute{\mathbf{V}}}_{\mathbf{k}}^{\mathbf{\sim}}\mathbf{,}\mathbf{B}_{\mathbf{k}}^{\mathbf{\sim}}\mathbf{\ \rbrack}\)
7. 　　　　\(\mathbf{if(k.eq.0)}\)
8. 　　　　　　\(\mathbf{{QR\ decomposition\ }}\mathbf{V}_{\mathbf{0}}^{\mathbf{\sim}}\mathbf{=}\mathbf{Q}_{\mathbf{0}}^{\mathbf{\sim}}\mathbf{R}_{\mathbf{0}}^{\mathbf{\sim}}\)
9. 　　　　　　\(\mathfrak{Q}_{\mathbf{0}}^{\mathbf{\sim}}\mathbf{=}\mathbf{Q}_{\mathbf{0}}^{\mathbf{\sim}}\)
10. 　　　　　　\(\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}}\)
11. 　　　　　　\(\mathcal{H}_{\mathbf{0}}\mathbf{=}\mathfrak{H}_{\mathbf{0}}^{\mathbf{\sim}}\)
12. 　　　　　　\(\mathbf{for\ o = 1,2\ldots,s}\)
13. 　　　　　　\(\mathbf{Givens\ rotation\ \lbrack o,}\mathcal{H}_{\mathbf{0}}\mathbf{,\zeta\rbrack}\)
14. 　　　　　　\(\mathbf{{if\ }}\left\| \mathbf{\zeta}_{\mathbf{o + 1}} \right\|\mathbf{\ /\ }\left\| \mathbf{b} \right\|\mathbf{{\ small\ enough\ then}}\)
15. 　　　　　　　　\(\mathbf{{quit}}\)
16. 　　　　　　\(\mathbf{{endif}}\)
17. 　　　　\(\mathbf{{else}}\)
18. 　　　　　　\({\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}}\)
19. 　　　　　　\({\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}}\)
20. 　　　　　　\(\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}}\)
21. 　　　　　　\(\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}\)
22. 　　　　　　\(\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}}\)
23. 　　　　　　\(\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}}\)
24. 　　　　　　\(\mathbf{h}_{\mathbf{k}}\mathbf{=}{\widetilde{\mathbf{\rho}}}_{\mathbf{k}}^{\mathbf{- 1}}\mathbf{\rho}_{\mathbf{k}}\mathbf{b}_{\mathbf{k}}\)
25. 　　　　　　\(\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}\)
26. 　　　　　　\(\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}\)
27. 　　　　　　\(\mathbf{for\ o = 1 + sk,\ldots,s(k + 1)}\)
28. 　　　　　　　　\(\mathbf{Givens\ rotation\ \lbrack o,}\mathcal{H}_{\mathbf{k}}\mathbf{,\zeta\rbrack}\)
29. 　　　　　　　　\(\mathbf{{if\ }}\left\| \mathbf{\zeta}_{\mathbf{o + 1}} \right\|\mathbf{\ /\ }\left\| \mathbf{b} \right\|\mathbf{{\ small\ enough\ then}}\)
30. 　　　　　　　　　　\(\mathbf{{quit}}\)
31. 　　　　　　　　\(\mathbf{{endif}}\)
32. 　　　　　　\(\mathbf{{end}}\)
33. 　　　　\(\mathbf{{endif}}\)
34. 　　\(\mathbf{{end}}\)
35. 　　\(\mathbf{Eigen\ Value\ problem\ \lbrack\mathfrak{H}_{t - 1}^{\sim},E,z\rbrack}\)
36. 　　\(\mathbf{Eigen Vectors \ X = \mathfrak{Q}_{t - 1}^{\sim}z}\)

2 前処理

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

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

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

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

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

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

　係数行列\(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\)の決定方法として,次の2とおりの条件のどちらを選択する.

• 条件1 対角成分が一致 \[A_{{ii}} = K_{{ii}}\ \ (i = 1,\ldots,n)\]

• 条件2 行内の要素和が一致 \[\sum_{j}^{}A_{{ij}} = \sum_{j}^{}K_{{ij}}\ \ (i = 1,\ldots,n)\]

対角行列\(D\)の具体的な計算方法は

• 条件1 対角成分が一致 \[D_{{ii}} = A_{{ii}} - \sum_{j < i}^{}{A_{{ij}}D_{jj}^{-1}A_{{ji}}}\ \ (i = 1,\ldots,n)\]

• 条件2 行内の要素和が一致 \[D_{{ii}} = A_{{ii}} - \sum_{j < i}^{}{\sum_{k > j}^{}{A_{{ij}}D_{{jj}}^{-1}A_{jk}}}\ \ (i = 1,\ldots,n)\]

となる.

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

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

2.5 加法的シュワルツ前処理 [1,2]

　前処理行列K として重なりのあるブロック対角行列を用いる.ブロック内で不完全LU 分解を用いる. 通信は必要となるが,計算は並列に行なわれる. PARCELでは加法的シュワルツ前処理の手法としてBASIC,RESTRICT,INTERPOLATE,NONEを選択できる.

BASIC

　\(r\)を重なりのある領域まで拡張させた\(Ks=r\)を解き,重なりのある領域の\(s\)を加算した結果を採用する.

RESTRICT

　\(r\)を重なりのある領域まで拡張させた\(Ks=r\)を解き,重なりのある領域の\(s\)を考慮しない結果を採用する.

INTERPOLATE

　\(r\)を重なりのある領域まで拡張させないで\(Ks=r\)を解き,重なりのある領域の\(s\)を加算した結果を採用する.

NONE

　\(r\)を重なりのある領域まで拡張させないで\(Ks=r\)を解き,重なりのある領域の\(s\)を考慮しない結果を採用する.

2.6 細分化前処理 [7]

　細分化前処理は不完全LU 分解でSIMD演算を適用可能となるように前処理行列を作成する方法である.ブロックヤコビ前処理及び加法シュワルツ前処理におけるブロック内で細分化対角ブロックを定義し、細分化対角ブロック内でのみ非対角要素を無視することでデータ依存性を排除してSIMD演算を実現する. 図6は左の9 \(\times\) 9のブロックにおいて3 \(\times\) 3の細分化対角ブロックを定義し、細分化対角ブロック内の非対角要素(白抜き)を無視することで3要素のベクトル処理をSIMD演算で行う.

ブロックヤコビ前処理	細分化ブロックヤコビ前処理
図6細分化前処理

2.7 ノイマン展開前処理 [2]

　係数行列\(A\)を対角スケーリングした行列\(\bar{A}\)を定義する.ノイマン展開前処理では前処理行列\(K\)として\(\bar{A}\)の逆行列をノイマン級数で近似する.

\[ \bar{A}^{-1}\simeq K^{-1}=[I+F+F^2+F^3+...] \]

ここで,\(\bar{A}=I-F\)となる.本手法は疎行列ベクトル積（SpMV）によって実装できるため,並列処理効率が高い.さらに,ブロックヤコビ前処理と組み合わせてSpMVのプロセス間通信を回避したブロックノイマン展開前処理も構築している.本前処理はGPU版のみで利用可能である.

3 QR分解

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

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

　グラムシュミットの直交化による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,2]

　誤差の低減を目的に古典グラムシュミット法のアルゴリズムを修正した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 Tall Skinny QR (TSQR)法 [3]

　TSQRは以下に示すハウスホルダーQRを基にした直交性の良いQR分解である. スレッド間はSequential TSQRで行い,プロセス間はParallel TSQRでQR分解を行う.

• アルゴリズム

1. \(\mathbf{for\ i = 1,\ldots,n\ do}\)
2. 　　\(\mathbf{y}_{\mathbf{i}}\mathbf{(i:m) = A(i:m,i) -}\left\| \mathbf{A(i:m,i)} \right\|\mathbf{e}_{\mathbf{i}}\mathbf{(i:m)}\)
3. 　　\(\mathbf{t}_{\mathbf{i}}\mathbf{=}\frac{\mathbf{2}}{\mathbf{<}\mathbf{y}_{\mathbf{i}}\left( \mathbf{i:m} \right)\mathbf{,\ }\mathbf{y}_{\mathbf{i}}\mathbf{(i:m) >}}\)
4. 　　\(\mathbf{Q}_{\mathbf{i}}\mathbf{=}\left( \mathbf{I -}\mathbf{t}_{\mathbf{i}}\mathbf{y}_{\mathbf{i}}\mathbf{(i:m)}{\mathbf{y}_{\mathbf{i}}\mathbf{(i:m)}}^{\mathbf{T}} \right)\)
5. 　　\(\mathbf{A(i:m,i) =}\mathbf{Q}_{\mathbf{i}}\mathbf{A(i:m,i)}\)
6. \(\mathbf{{enddo}}\)

Sequential TSQR		Parallel TSQR

3.4 コレスキーQR法 [8]

　コレスキー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}}\)

3.5 コレスキーQR2法 [9]

　コレスキーQRにより得られた直交行列に対して再度コレスキーQRを行うアルゴリズムである. コレスキーQRを2回実行することで直交性を改善させる.

• アルゴリズム

1. 　\(\mathbf{B}_{\mathbf{1}}\mathbf{=}\mathbf{V}^{\mathbf{T}}\mathbf{V}\)
2. 　\(\mathbf{R}_{\mathbf{1}}^{\mathbf{T}}\mathbf{R}_{\mathbf{1}}\mathbf{= Cholesky\ decomposition(}\mathbf{B}_{\mathbf{1}}\mathbf{)}\)
3. 　\(\mathbf{Q}_{\mathbf{1}}\mathbf{= V}\mathbf{R}_{\mathbf{1}}^{\mathbf{- 1}}\)
4. 　\(\mathbf{B}_{\mathbf{2}}\mathbf{=}\mathbf{Q}_{\mathbf{1}}^{\mathbf{T}}\mathbf{Q}_{\mathbf{1}}\)
5. 　\(\mathbf{R}_{\mathbf{2}}^{\mathbf{T}}\mathbf{R}_{\mathbf{2}}\mathbf{= Cholesky\ decomposition(}\mathbf{B}_{\mathbf{2}}\mathbf{)}\)
6. 　\(\mathbf{Q}_{\mathbf{2}}\mathbf{=}\mathbf{Q}_{\mathbf{1}}\mathbf{R}_{\mathbf{2}}^{\mathbf{- 1}}\)
7. 　\(\mathbf{V =}\mathbf{Q}_{\mathbf{1}}\mathbf{R}_{\mathbf{1}}\mathbf{=}\mathbf{Q}_{\mathbf{2}}\mathbf{R}_{\mathbf{2}}\mathbf{R}_{\mathbf{1}}\mathbf{= QR}\)
8. 　\(\mathbf{Q =}\mathbf{Q}_{\mathbf{2}}\)
9. 　\(\mathbf{R =}\mathbf{R}_{\mathbf{2}}\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\}\)

4.2 Diagonal (DIA)形式

　DIA形式(圧縮対角形式)は,対角行単位に格納する方法である.対角行の値と対角からのオフセットを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プロセス

対角要素　　　　::　\(diaA=\{0,0,e,h,0,0,0,i,a,d,f,j,b,0,g,0,c,0,0,0\}\)

オフセット　　　::　\(offset=\{-2,-1,0,1,2\}\)

対角要素の本数　::　\(nnd=5\)

2プロセスの場合

• ランク0プロセス

対角要素　　　　::　\(diaA=\{a,d,b,0,c,0\}\)

オフセット　　　::　\(offset=\{0,1,2\}\)

対角要素の本数　::　\(nnd=3\)

• ランク1プロセス

対角要素　　　　::　\(diaA=\{e,h,0,i,f,j,g,0\}\)

オフセット　　　::　\(offset=\{-2,-1,0,1\}\)

対角要素の本数　::　\(nnd=4\)

4.3 Domain Decomposition Method (DDM)形式

4.3.1 多次元領域分割におけるパラメータ

　DDM形式(多次元領域分割形式)は,多次元領域分割した差分法等のステンシル計算に基づくシミュレーションを想定してDIA形式を拡張した行列形式である. 通信するプロセスのリストをシミュレーション側の多次元領域分割と一致するように指定することにより,1次元,2次元,ないしは,3次元の構造格子で与えられる行列の並列処理に必要となる通信を最小化することができる. 関係するパラメータ及び指定法を以下に説明する.

• 変数：num_neighbor_ranks

　多次元領域分割における隣接プロセスの数を意味する. \(m\)をプロセスの分割軸の数として,3\(^{m-1}\)で与えられる. 例えば,2軸分割の場合は8,3軸分割の場合は26 である.

• 1次元配列：neighbor_ranks ( num_neighbor_ranks )

　MPIコミュニケータ（1次元）上での隣接するプロセス番号のリストを指定する. なお,計算領域の境界で隣接プロセスが存在しない方向については負のインデックスを指定する. PARCELにおけるneighbor_ranksの指定法の例を以下に示す.ここで,npe_x, npe_y, npe_z は,それぞれx,y,z 方向の構造格子の分割数,rank_x, rank_y, rank_z はプロセス番号の各方向におけるランクである.ただし,1次元化したインデックスiirankの定義はシミュレーションとそろえること.

• 2次元配列：ioff_grids (nnd, ncomp_grids )

　差分法等のステンシル計算が参照する位置を構造格子における各方向に指定する. 例えば,3次元ポアソン方程式の3点差分からなる7点ステンシルであれば,差分計算で使用する東西南北上下の格子点に対する基準点からのオフセットを以下のように指定する.

• 行列構造

1次元ポアソン方程式の3点差分の場合.

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

1プロセスの場合

• ランク0プロセス

対角要素　　　　::　\(val\_dia=\{a,d,g,j,b,e,h,0,0,c,f,i\}\)

オフセット　　　::　\(ioff\_dia=\{0,0,-1\}\)

対角要素の本数　::　\(nnd=3\)

ステンシル成分　::

2プロセスの場合

• ランク0プロセス

対角要素　　　　::　\(val\_dia=\{a,d,b,e,0,c\}\)

オフセット　　　::　\(ioff\_dia=\{0,0,-1\}\)

対角要素の本数　::　\(nnd=3\)

ステンシル成分　::

• ランク1プロセス

対角要素　　　　::　\(val\_dia=\{g,j,h,0,f,i\}\)

オフセット　　　::　\(ioff\_dia=\{0,0,0\}\)

対角要素の本数　::　\(nnd=3\)

ステンシル成分　::

5 インストール

　PARCELの導入方法について説明する.

5.1 コンパイル方法

　PARCELのコンパイルにはC/C++コンパイラ及びプリプロセッサ対応のFortranコンパイラ,MPI,LAPACK,OpenMPライブラリを必要とする.GPU版を利用する場合はCUDAが追加で必要となる. gfortranでコンパイルする場合はコンパイルオプションに-cppを追加することでプリプロセッサが適用される. コンパイルオプションはarch/make_configを編集することで変更する. make_configの記述例を以下に示す. コンパイル方法は展開したアーカイブのトップディレクトリでmakeコマンドを実行する. 続いて”make install”を実行することでINSTALLDIRに指定したディレクトリにexample,include,libディレクトリが作成される. libディレクトにparcel_sparse.aが作成されていればコンパイル成功である.exampleディレクトリでmakeコマンドを実行することでサンプルプログラムがコンパイルされる. 20億次元以上の大規模問題では64ビット整数が必要となるため64ビット整数型でコンパイルする必要がある. GPU版PARCELを利用する場合はCUDAFLAGをONにしてコンパイルする.CUDAAWAREをONにすることでCUDA AWARE MPI環境に対応するライブラリが作成される.GPU版PARCELではcuBLAS,cuSPARSEを利用するためライブラリをリンクする必要がある.

• make_configの記述例

  CUDAFLAG  = ON
  CUDAAWARE = ON
  CC     = mpiicc
  FC     = mpiifort
 
  CFLAGS  = -O3
  FFLAGS  = -O3 -fpp -qopenmp -xHost  -mcmodel=large
 
  FP_MODE = -fp-model precise
 
  INCLUDEDIR =
  LD_MPI     = -lmpifort
  LD_FORT    = -lifcore
  LD_LAPACK  = -mkl
  LD_CUDA    = -lcudart -lcublas -lcusparse
  NVCFLAGS  = \
            -arch=sm_80 -ccbin=${CC} -O3 -restrict  \
            --expt-extended-lambda -use_fast_math   \
            --default-stream per-thread  \
            --nvlink-options -Werror     \
            --maxrregcount=128  \
            --extra-device-vectorization  \
            --restrict \
             "-std=c++11"
 
  #FDEFINE = -DPARCEL_INT8
  INSTALLDIR = ../PARCEL_1.2

5.2 使用方法

　libディレクトリに生成されるparcel_sparse.aファイルをリンクすることでPARCELのサブルーチンを使用可能となる.

6 PARCELサブルーチン

　PARCELに整備されているサブルーチンの呼び出しインターフェースについて説明する. 以下の実装は各ルーチンで共通に使用されている.

• 演算精度

　parcel_dから始まる関数名は倍精度、parcel_ddから始まる関数名は4倍精度（Double-Double形式）で実装されている. GPU版は4倍精度に未対応である.

• GPU版

　GPU版はFortranのtype(c_ptr)にGPU側のメモリを割り当てている. プロセスが利用するGPUを指定するためにクリロフ部分空間法のサブルーチンを実行する前にサブルーチンPARCEL_Initを呼ぶ必要がある. 入力する係数行列A,右辺ベクトルb,初期解ベクトルxについてもGPU側のメモリを割り当てたtype(c_ptr)を利用する. CRS形式の疎行列ベクトル積にはcuSPARSE(cusparseSpMV)を利用している. GPU版は細分化前処理に未対応であるがGPUに適したノイマン展開前処理が利用可能である.

• 通信隠蔽

　CPU版では通信隠蔽について3種類の実装が利用可能である. iovlflag=0では通信隠蔽を行わない. iovlflag=1では全ての非同期通信処理を同時発行した後に計算処理を実行し,通信と演算を同時に処理する. iovlflag=0では各方向の非同期通信と袖領域の計算処理を同期を取りながら順番に隠蔽処理することで通信バッファを節約する. GPU版は通信隠蔽に未対応である.

• ブロックヤコビ前処理,加法シュワルツ前処理のブロック数

　CPU版ではprecon_thblockに0以下もしくはプロセスあたりのスレッド数を上回る値を設定するとデフォルト値としてプロセスあたりのスレッド数が設定される. スレッド数を減らしてブロックサイズを拡大することで前処理の精度が向上するが並列処理性能は低下する. GPU版では1プロセス1ブロックをcuSPARSEのILU分解(cusparseDcsrsv2_solve)により計算する.

• 細分化前処理で独立に計算する対角ブロックのサイズ

　independ_nvecに0以下の値を設定した場合には細分化処理は適用されない. A64FXでは、富士通コンパイラのSIMD処理とソフトウェアパイプライニングに必要なサイズとしてindepend_nvec=300程度が推奨される. independ_nvecを大きくすると処理性能が向上するが前処理の精度は低下する. independ_nvecが各プロセスが担当するベクトルの大きさnを超えた場合にはポイントヤコビ前処理と等価な処理となる.

• SpMVのCPUスレッド最適化

　SpMVのスレッド並列処理はスレッド間のキャッシュ再利用を意識したサイクリックループ分割を適用している. この際のチャンクサイズをnBlockで指定する. nBlockに0以下の値を設定した場合にはデフォルト値のnBlock=2000が適用される.

• 通信削減型アーノルディ法における固有ベクトルの格納方法

　i番目の固有値が実数であれば対応する固有ベクトルはEvec(n*(i-1)+1:n*i)に格納される. i番目の固有値が複素数であり,i+1番目の固有値がi番目の固有値の共役複素数の場合のi,i+1番目の固有ベクトルは次のように格納される. GPU版は固有値計算に未対応である.

6.1 PARCEL_Init

• 呼び出し方法

　call PARCEL_Init(myrank)

• 機能

　PARCELの初期化処理を行う.

6.2 PARCEL_Free

• 呼び出し方法

　call PARCEL_Free()

• 機能

　PARCELの終了処理を行う.

6.3 parcel_dcg (CPU/GPU)

• 呼び出し方法

call parcel_dcg( icomm, vecx,vecb, n,gn,nnz,istart, crsA,crsRow_ptr,crsCol, ipreflag,ilu_method,addL,iflagAS, itrmax,rtolmax, reshistory, iovlflag, precon_thblock, Neuman_itrmax, independ_nvec, nBlock, iret )

• 機能

　共役勾配法(CG法)により,連立一次方程式Ax=bを解く.非零要素の値をCRS形式で格納.

6.4 parcel_dbicgstab (CPU/GPU)

• 呼び出し方法

call parcel_dbicgstab( icomm, vecx,vecb, n,gn,nnz,istart, crsA,crsRow_ptr,crsCol, ipreflag,ilu_method,addL,iflagAS, itrmax,rtolmax, reshistory, iovlflag, precon_thblock, Neuman_itrmax, independ_nvec, nBlock, iret )

• 機能

　安定化双共役勾配法(Bi-CGSTAB法)により,連立一次方程式Ax=bを解く.非零要素の値をCRS形式で格納.

6.5 parcel_dgmres (CPU/GPU)

• 呼び出し方法

call parcel_dgmres( icomm, vecx,vecb, n,gn,nnz,istart, crsA,crsRow_ptr,crsCol, ipreflag,ilu_method,addL,iflagAS, itrmax,rtolmax, reshistory, iovlflag, precon_thblock, Neuman_itrmax, independ_nvec, nBlock, gmres_m,gmres_GSflag, iret )

• 機能

　一般化最小残差法(GMRES(m)法)により,連立一次方程式Ax=bを解く.非零要素の値をCRS形式で格納.

6.6 parcel_dcagmres (CPU/GPU)

• 呼び出し方法

call parcel_dcagmres( icomm, vecx,vecb, n,gn,nnz,istart, crsA,crsRow_ptr,crsCol, ipreflag,ilu_method,addL,iflagAS, itrmax,rtolmax, reshistory, iovlflag, precon_thblock, Neuman_itrmax, independ_nvec, nBlock, cagmres_sstep, cagmres_tstep, cagmres_basis, cagmres_QRflag, iret )

• 機能

　通信削減型一般化最小残差法(CA-GMRES(s,t)法)により,連立一次方程式Ax=bを解く.非零要素の値をCRS形式で格納.

6.7 parcel_dcbcg (CPU/GPU)

• 呼び出し方法

call parcel_dcbcg( icomm, vecx,vecb, n,gn,nnz,istart, crsA,crsRow_ptr,crsCol, ipreflag,ilu_method,addL,iflagAS, itrmax,rtolmax, reshistory, iovlflag, precon_thblock, Neuman_itrmax, independ_nvec, nBlock, cbcg_kstep,cbcg_Eigenflag,power_method_itrmax, caarnoldi_sstep,caarnoldi_tstep, caarnoldi_basis,caarnoldi_QRflag, iret )

• 機能

　チェビシェフ基底共役勾配法(CBCG法)により,連立一次方程式Ax=bを解く.非零要素の値をCRS形式で格納.

6.8 parcel_dcaarnoldi (CPU)

• 呼び出し方法

call parcel_dcaarnoldi( icomm, vecx, vecb, n, gn, nnz, istart, crsA, crsRow_ptr, crsCol, ipreflag, ilu_method, addL, iflagAS, itrmax, iovlflag, precon_thblock, independ_nvec, nBlock, caarnoldi_sstep, caarnoldi_tstep, caarnoldi_basis, caarnoldi_QRflag, Evalr, Evali, Evec, Eerr, EmaxID, EminID, iret )

• 機能

　通信削減型アーノルディ法(CA-Arnoldi((s,t)法)により,固有値問題を解く.非零要素の値をCRS形式で格納.

6.9 parcel_dcg_dia (CPU/GPU)

• 呼び出し方法

call parcel_dcg_dia( icomm, vecx,vecb, n,gn,istart, diaA,offset,nnd, ipreflag,ilu_method,addL,iflagAS, itrmax,rtolmax, reshistory, iovlflag, precon_thblock, Neuman_itrmax, independ_nvec, nBlock, iret )

• 機能

　共役勾配法(CG法)により,連立一次方程式Ax=bを解く.非零要素の値をDIA形式で格納.

6.10 parcel_dbicgstab_dia (CPU/GPU)

• 呼び出し方法

call parcel_dbicgstab_dia( icomm, vecx,vecb, n,gn,istart, diaA,offset,nnd, ipreflag,ilu_method,addL,iflagAS, itrmax,rtolmax, reshistory, iovlflag, precon_thblock, Neuman_itrmax, independ_nvec, nBlock, iret )

• 機能

　安定化双共役勾配法(Bi-CGSTAB法)により,連立一次方程式Ax=bを解く.非零要素の値をDIA形式で格納.

6.11 parcel_dgmres_dia (CPU/GPU)

• 呼び出し方法

call parcel_dgmres_dia( icomm, vecx,vecb, n,gn,istart, diaA,offset,nnd, ipreflag,ilu_method,addL,iflagAS, itrmax,rtolmax, reshistory, iovlflag, precon_thblock, Neuman_itrmax, independ_nvec, nBlock, gmres_m,gmres_GSflag, iret )

• 機能

　一般化最小残差法(GMRES(m)法)により,連立一次方程式Ax=bを解く.非零要素の値をDIA形式で格納.

6.12 parcel_dcagmres_dia (CPU/GPU)

• 呼び出し方法

call parcel_dcagmres_dia( icomm, vecx,vecb, n,gn,istart, diaA,offset,nnd, ipreflag,ilu_method,addL,iflagAS, itrmax,rtolmax, reshistory, iovlflag, precon_thblock, Neuman_itrmax, independ_nvec, nBlock, cagmres_sstep, cagmres_tstep, cagmres_basis, cagmres_QRflag, iret )

• 機能

　通信削減型一般化最小残差法(CA-GMRES(s,t)法)により,連立一次方程式Ax=bを解く.非零要素の値をDIA形式で格納.

6.13 parcel_dcbcg_dia (CPU/GPU)

• 呼び出し方法

call parcel_dcbcg_dia( icomm, vecx,vecb, n,gn,istart, diaA,offset,nnd, ipreflag,ilu_method,addL,iflagAS, itrmax,rtolmax, reshistory, iovlflag, precon_thblock, Neuman_itrmax, independ_nvec, nBlock, cbcg_kstep,cbcg_Eigenflag,power_method_itrmax, caarnoldi_sstep,caarnoldi_tstep, caarnoldi_basis,caarnoldi_QRflag, iret )

• 機能

　チェビシェフ基底共役勾配法(CBCG法)により,連立一次方程式Ax=bを解く.非零要素の値をDIA形式で格納.

6.14 parcel_dcaarnoldi_dia (CPU)

• 呼び出し方法

call parcel_dcaarnoldi_dia( icomm, vecx, vecb, n, gn, istart, diaA,offset,nnd, ipreflag, ilu_method, addL, iflagAS, itrmax, iovlflag, precon_thblock, independ_nvec, nBlock, caarnoldi_sstep, caarnoldi_tstep, caarnoldi_basis, caarnoldi_QRflag, Evalr, Evali, Evec, Eerr, EmaxID, EminID, iret )

• 機能

　通信削減型アーノルディ法(CA-Arnoldi((s,t)法)により,固有値問題を解く.非零要素の値をCRS形式で格納.

6.15 parcel_dcg_ddm (CPU/GPU)

• 呼び出し方法

call parcel_dcg_ddm( icomm, x,b, m, nnd,ioff_dia,val_dia, num_neighbor_ranks,neighbor_ranks,margin,ncomp_grids, ndiv_grids,num_grids,ioff_grids, div_direc_th, ipreflag, addL,iflagAS, itrmax,rtolmax,reshistory, iovlflag, precon_thblock, Neuman_itrmax, independ_nvec, nBlock, iret )

• 機能

　共役勾配法(CG法)により,連立一次方程式Ax=bを解く.非零要素の値をDDM形式で格納.

6.16 parcel_dbicgstab_ddm (CPU/GPU)

• 呼び出し方法

call parcel_dbicgstab_ddm( icomm, x,b, m, nnd,ioff_dia,val_dia, num_neighbor_ranks,neighbor_ranks,margin,ncomp_grids, ndiv_grids,num_grids,ioff_grids, div_direc_th, ipreflag, addL,iflagAS, itrmax,rtolmax,reshistory, iovlflag, precon_thblock, Neuman_itrmax, independ_nvec, nBlock, iret )

• 機能

　安定化双共役勾配法(Bi-CGSTAB法)により,連立一次方程式Ax=bを解く.非零要素の値をDDM形式で格納.

6.17 parcel_dgmres_ddm (CPU/GPU)

• 呼び出し方法

call parcel_dgmres_ddm( icomm, x,b, m, nnd,ioff_dia,val_dia, num_neighbor_ranks,neighbor_ranks,margin,ncomp_grids, ndiv_grids,num_grids,ioff_grids, div_direc_th, ipreflag, addL,iflagAS, itrmax,rtolmax,reshistory, iovlflag, precon_thblock, Neuman_itrmax, independ_nvec, nBlock, gmres_m, gmres_GSflag, iret )

• 機能

　一般化最小残差法(GMRES(m)法)により,連立一次方程式Ax=bを解く.非零要素の値をDDM形式で格納.

6.18 parcel_dcagmres_ddm (CPU/GPU)

• 呼び出し方法

call parcel_dcagmres_ddm( icomm, x,b, m, nnd,ioff_dia,val_dia, num_neighbor_ranks,neighbor_ranks,margin,ncomp_grids, ndiv_grids,num_grids,ioff_grids, div_direc_th, ipreflag, addL,iflagAS, itrmax,rtolmax,reshistory, iovlflag, precon_thblock, Neuman_itrmax, independ_nvec, nBlock, cagmres_sstep, cagmres_tstep, cagmres_basis, cagmres_QRflag, iret )

• 機能

　通信削減型一般化最小残差法(CA-GMRES(s,t)法)により,連立一次方程式Ax=bを解く.非零要素の値をDDM形式で格納.

6.19 parcel_dcbcg_ddm (CPU/GPU)

• 呼び出し方法

call parcel_dcbcg_ddm( icomm, x,b, m, nnd,ioff_dia,val_dia, num_neighbor_ranks,neighbor_ranks,margin,ncomp_grids, ndiv_grids,num_grids,ioff_grids, div_direc_th, ipreflag, addL,iflagAS, itrmax,rtolmax,reshistory, iovlflag, precon_thblock, Neuman_itrmax, independ_nvec, nBlock, cbcg_kstep, cbcg_Eigenflag, power_method_itrmax, caarnoldi_sstep, caarnoldi_tstep, caarnoldi_basis, caarnoldi_QRflag, iret )

• 機能

　チェビシェフ基底共役勾配法(CBCG法)により,連立一次方程式Ax=bを解く.非零要素の値をDDM形式で格納.

6.20 parcel_dcaarnoldi_ddm (CPU)

• 呼び出し方法

call parcel_dcaarnoldi_ddm( icomm, x,b, m, nnd,ioff_dia,val_dia, num_neighbor_ranks,neighbor_ranks,margin,ncomp_grids, ndiv_grids,num_grids,ioff_grids, div_direc_th, ipreflag, addL,iflagAS, iovlflag, precon_thblock, independ_nvec, nBlock, caarnoldi_sstep, caarnoldi_tstep, caarnoldi_basis, caarnoldi_QRflag, Evalr, Evali, Evec, Eerr, EmaxID, EminID, iret )

• 機能

　通信削減型アーノルディ法(CA-Arnoldi((s,t)法)により,連立一次方程式Ax=bを解く.非零要素の値をDDM形式で格納.

6.21 parcel_ddcg (CPU)

• 呼び出し方法

call parcel_ddcg( icomm, vecx_hi,vecx_lo, vecb_hi,vecb_lo, n,gn,nnz,istart, crsA_hi,crsA_lo,crsRow_ptr,crsCol, ipreflag,ilu_method,addL,iflagAS, itrmax,rtolmax, reshistory, iovlflag, precon_thblock, independ_nvec, nBlock, precision_A, precision_b, precision_x, precision_precon, iret )

• 機能

　共役勾配法(CG法)により,連立一次方程式Ax=bを解く.非零要素の値をCRS形式で格納.

6.22 parcel_ddbicgstab (CPU)

• 呼び出し方法

call parcel_ddbicgstab( icomm, vecx_hi,vecx_lo, vecb_hi,vecb_lo, n,gn,nnz,istart, crsA_hi,crsA_lo,crsRow_ptr,crsCol, ipreflag,ilu_method,addL,iflagAS, itrmax,rtolmax, reshistory, iovlflag, precon_thblock, independ_nvec, nBlock, precision_A, precision_b, precision_x, precision_precon, iret )

• 機能

　安定化双共役勾配法(Bi-CGSTAB法)により,連立一次方程式Ax=bを解く.非零要素の値をCRS形式で格納.

6.23 parcel_ddgmres (CPU)

• 呼び出し方法

call parcel_ddgmres( icomm, vecx_hi,vecx_lo, vecb_hi,vecb_lo, n,gn,nnz,istart, crsA_hi,crsA_lo,crsRow_ptr,crsCol, ipreflag,ilu_method,addL,iflagAS, itrmax,rtolmax, reshistory, iovlflag, precon_thblock, independ_nvec, nBlock, gmres_m,gmres_GSflag, precision_A, precision_b, precision_x, precision_precon, iret )

• 機能

　一般化最小残差法(GMRES(m)法)により,連立一次方程式Ax=bを解く.非零要素の値をCRS形式で格納.

6.24 parcel_ddcg_dia (CPU)

• 呼び出し方法

call parcel_ddcg_dia( icomm, vecx_hi,vecx_lo, vecb_hi,vecb_lo, n,gn,istart, diaA_hi,diaA_lo,offset,nnd, ipreflag,ilu_method,addL,iflagAS, itrmax,rtolmax, reshistory, iovlflag, precon_thblock, independ_nvec, nBlock, precision_A, precision_b, precision_x, precision_precon, iret )

• 機能

　共役勾配法(CG法)により,連立一次方程式Ax=bを解く.非零要素の値をDIA形式で格納.

6.25 parcel_ddbicgstab_dia (CPU)

• 呼び出し方法

call parcel_ddbicgstab_dia( icomm, vecx_hi,vecx_lo, vecb_hi,vecb_lo, n,gn,istart, diaA_hi,diaA_lo,offset,nnd, ipreflag,ilu_method,addL,iflagAS, itrmax,rtolmax, reshistory, iovlflag, precon_thblock, independ_nvec, nBlock, precision_A, precision_b, precision_x, precision_precon, iret )

• 機能

　安定化双共役勾配法(Bi-CGSTAB法)により,連立一次方程式Ax=bを解く.非零要素の値をDIA形式で格納.

6.26 parcel_ddgmres_dia (CPU)

• 呼び出し方法

call parcel_ddgmres_dia( icomm, vecx_hi,vecx_lo, vecb_hi,vecb_lo, n,gn,istart, diaA_hi,diaA_lo,offset,nnd, ipreflag,ilu_method,addL,iflagAS, itrmax,rtolmax, reshistory, iovlflag, precon_thblock, independ_nvec, nBlock, precision_A, precision_b, precision_x, precision_precon, iret )

• 機能

　一般化最小残差法(GMRES(m)法)により,連立一次方程式Ax=bを解く.非零要素の値をDIA形式で格納.

6.27 parcel_ddcg_ddm (CPU)

• 呼び出し方法

call parcel_ddcg_ddm( icomm, x_hi,x_lo, b_hi,b_lo, m, nnd,ioff_dia, val_dia_hi,val_dia_lo, num_neighbor_ranks,neighbor_ranks,margin,ncomp_grids, ndiv_grids,num_grids,ioff_grids, div_direc_th, ipreflag, addL,iflagAS, itrmax,rtolmax,reshistory, iovlflag, precon_thblock, independ_nvec, nBlock, precision_A, precision_b, precision_x, precision_precon, iret )

• 機能

　共役勾配法(CG法)により,連立一次方程式Ax=bを解く.非零要素の値をDDM形式で格納.

6.28 parcel_ddbicgstab_ddm (CPU)

• 呼び出し方法

call parcel_ddbicgstab_ddm( icomm, x_hi,x_lo, b_hi,b_lo, m, nnd,ioff_dia, val_dia_hi,val_dia_lo, num_neighbor_ranks,neighbor_ranks,margin,ncomp_grids, ndiv_grids,num_grids,ioff_grids, div_direc_th, ipreflag, addL,iflagAS, itrmax,rtolmax,reshistory, iovlflag, precon_thblock, independ_nvec, nBlock, precision_A, precision_b, precision_x, precision_precon, iret )

• 機能

　安定化双共役勾配法(Bi-CGSTAB法)により,連立一次方程式Ax=bを解く.非零要素の値をDDM形式で格納.

6.29 parcel_ddgmres_ddm (CPU)

• 呼び出し方法

call parcel_ddgmres_ddm( icomm, x_hi,x_lo, b_hi,b_lo, m, nnd,ioff_dia, val_dia_hi,val_dia_lo, num_neighbor_ranks,neighbor_ranks,margin,ncomp_grids, ndiv_grids,num_grids,ioff_grids, div_direc_th, ipreflag, addL,iflagAS, itrmax,rtolmax,reshistory, iovlflag, precon_thblock, independ_nvec, nBlock, gmres_m, gmres_GSflag, precision_A, precision_b, precision_x, precision_precon, iret )

• 機能

　一般化最小残差法(GMRES(m)法)により,連立一次方程式Ax=bを解く.非零要素の値をDDM形式で格納.

6.30 parcel_dqr (CPU)

• 呼び出し方法

call parcel_dqr(n,s,icomm,V,Q,R,iQRflag)

• 機能

　m行s列の行列をQR分解する.ここで,MPIコミュニケータicommに属するMPIプロセスに対するnの総和がmとなる.

7 Fortran用CUDAサポートサブルーチン

　Fortran用に整備されているCUDAサポートサブルーチンの呼び出しインターフェースについて説明する. 本機能はCUDA FortranやOpenACCを用いて代替することも可能である.

7.1 cudaMalloc_FP64 (GPU)

• 呼び出し方法

call cudaMalloc_FP64(a_d,n)

• 機能

　double precision型のGPUメモリをtype(devicemem)構造体のtype(c_ptr)メンバmemに割り当てる.

7.2 cudaMalloc_PARCEL_INT (GPU)

• 呼び出し方法

call cudaMalloc_PARCEL_INT(a_d,n)

• 機能

　32ビット整数バージョンではinteger*4型,64ビット整数バージョンではinteger*8型のGPUメモリをtype(devicemem)構造体のtype(c_ptr)メンバmemに割り当てる.

7.3 cudaMalloc_INT (GPU)

• 呼び出し方法

call cudaMalloc_INT(a_d,n)

• 機能

　integer*4型のGPUメモリをtype(devicemem)構造体のtype(c_ptr)メンバmemに割り当てる.

7.4 cudaFree (GPU)

• 呼び出し方法

call cudaFree(a_d)

• 機能

　GPUメモリを解放する.

7.5 cudaMemcpy_h2d_FP64 (GPU)

• 呼び出し方法

call cudaMemcpy_h2d_FP64(a_h,a_d,n,offset_h,offset_d)

• 機能

　FP64データをCPUメモリからGPUメモリにコピー.

7.6 cudaMemcpy_h2d_PARCEL_INT (GPU)

• 呼び出し方法

call cudaMemcpy_h2d_PARCEL_INT(a_h,a_d,n,offset_h,offset_d)

• 機能

　32ビット整数バージョンではinteger*4型,64ビット整数バージョンではinteger*8型のデータをCPUメモリからGPUメモリにコピー.

7.7 cudaMemcpy_h2d_INT (GPU)

• 呼び出し方法

call cudaMemcpy_h2d_INT(a_h,a_d,n,offset_h,offset_d)

• 機能

　integer*4データをCPUメモリからGPUメモリにコピー.

7.8 cudaMemcpy_d2h_FP64 (GPU)

• 呼び出し方法

call cudaMemcpy_d2h_FP64(a_d,a_h,n,offset_d,offset_h)

• 機能

　FP64データをGPUメモリからCPUメモリにコピー.

7.9 cudaMemcpy_d2h_PARCEL_INT (GPU)

• 呼び出し方法

call cudaMemcpy_d2h_PARCEL_INT(a_d,a_h,n,offset_d,offset_h)

• 機能

　32ビット整数バージョンではinteger*4型,64ビット整数バージョンではinteger*8型のデータをGPUメモリからCPUメモリにコピー.

7.10 cudaMemcpy_d2h_INT (GPU)

• 呼び出し方法

call cudaMemcpy_d2h_INT(a_d,a_h,n,offset_d,offset_h)

• 機能

　integer*4データをGPUメモリからCPUメモリにコピー.

8 サブルーチン使用方法(Fortran)

　PARCELのサブルーチン使用方法についてCG法のFortranサンプルコードを例にして説明する. ここで,プリプロセッサでCUDA_SOLVERを有効にすればCUDA版のサンプルコードになる.

8.1 CRS形式 (Fortran)

　CRS形式のサンプルコードを以下に示す. make_matrix_CRSはPARCELの行列分割に対応したCRS形式の行列を生成する任意のサブルーチンとする.

  program main
    use mpi
    use PARCEL_integer
    use krylov
    use PARCEL_wrapper
    implicit none
    integer n,gn,nnz,istart
    real*8,allocatable  :: crsA(:)
    integer,allocatable :: crsRow_ptr(:),crsCol(:)
    real*8,allocatable  :: vecx(:)
    real*8,allocatable  :: vecb(:)
    integer itrmax
    real*8  rtolmax
    real*8, allocatable :: reshistory(:)
    integer ipreflag
    integer ILU_method
    integer iflagAS
    integer addL
    integer iovlflag
    integer precon_thblock
    integer independ_nvec
    integer nblock
    integer iret
    integer myrank
    integer ierr
 #ifdef CUDA_SOLVER
    type(devicemem) cu_vecx
    type(devicemem) cu_vecb
    type(devicemem) cu_crsA
    type(devicemem) cu_crsRow_ptr
    type(devicemem) cu_crsCol
    integer(kind=intType),parameter :: offset_zero=0
 #endif
    call MPI_Init(ierr)
    call MPI_Comm_rank(MPI_COMM_WORLD,myrank,ierr)
    call PARCEL_Init(myrank)
    call make_matrix_CRS(n,gn,nnz,istart,crsA,crsRow_ptr,crsCol)
    allocate(vecx(n))
    allocate(vecb(n))
    allocate(reshistory(itrmax))
    ipreflag=0
    ILU_method=1
    addL=0
    iflagAS=1
    itrmax=100
    rtolmax=1.0d-8
    iovlflag=0
    precon_thblock=-1
    independ_nvec=-1
    nblock = 2000
    vecb=1.0d0
    vecx=1.0d0
 #ifdef CUDA_SOLVER
    call cudaMalloc_FP64(cu_vecx,n)
    call cudaMalloc_FP64(cu_vecb,n)
    call cudaMalloc_FP64(cu_crsA,nnz)
    call cudaMalloc_PARCEL_INT(cu_crsRow_ptr,n+1)
    call cudaMalloc_PARCEL_INT(cu_crsCol,nnz)
    call cudaMemcpy_h2d_fp64(vecx,cu_vecx,n,offset_zero,offset_zero)
    call cudaMemcpy_h2d_fp64(vecb,cu_vecb,n,offset_zero,offset_zero)
    call cudaMemcpy_h2d_fp64(crsMat%crsA,cu_crsA,nnz,offset_zero,offset_zero)
    call cudaMemcpy_h2d_PARCEL_INT(crsMat%crsRow_ptr,cu_crsRow_ptr,n+1,offset_zero,offset_zero)
    call cudaMemcpy_h2d_PARCEL_INT(crsMat%crsCol,cu_crsCol,nnz,offset_zero,offset_zero)
 #endif
    call parcel_dcg( &
         MPI_COMM_WORLD, &
 #ifdef CUDA_SOLVER
         cu_vecx%mem,cu_vecb%mem, &
         n,gn,nnz,istart, &
         cu_crsA%mem,cu_crsRow_ptr%mem,cu_crsCol%mem,&
 #else
         vecx,vecb, &
         n,gn,nnz,istart, &
         crsA,crsRow_ptr,crsCol, &
 #endif
         ipreflag,ILU_method,addL,iflagAS, &
         itrmax,rtolmax, &
         reshistory, &
         iovlflag,iret &
         precon_thblock, independ_nvec, &
         nblock, &
         iret &
         )
    call PARCEL_Free()
    call MPI_Finalize(ierr)
 #ifdef CUDA_SOLVER
    call cudaMemcpy_d2h_fp64(cu_vecx,vecx,n,offset_zero,offset_zero)
    call cudaFree(cu_vecx)
    call cudaFree(cu_vecb)
    call cudaFree(cu_crsA)
    call cudaFree(cu_crsRow_ptr)
    call cudaFree(cu_crsCol)
 #endif
    deallocate(vecx)
    deallocate(vecb)
    deallocate(reshistory)
  end program main

8.2 DIA形式 (Fortran)

　DIA形式のサンプルコードを以下に示す. make_matrix_DIAはPARCELのブロック分割に対応したDIA形式の行列を生成する任意のサブルーチンとする.

  program main
    use PARCEL_integer
    use krylov
    use PARCEL_wrapper
    implicit none
    integer n,gn,nnd,istart
    real*8,allocatable  :: diaA(:)
    integer,allocatable :: offset(:)
    real*8,allocatable  :: vecx(:)
    real*8,allocatable  :: vecb(:)
    integer itrmax
    real*8  rtolmax
    real*8, allocatable :: reshistory(:)
    integer ipreflag
    integer ILU_method
    integer iflagAS
    integer addL
    integer iovlflag
    integer precon_thblock
    integer independ_nvec
    integer nblock
    integer iret
    integer myrank
    integer ierr
 #ifdef CUDA_SOLVER
    type(devicemem) cu_vecx
    type(devicemem) cu_vecb
    type(devicemem) cu_diaA
    type(devicemem) cu_offset
    integer(kind=intType),parameter :: offset_zero=0
 #endif
    call MPI_Init(ierr)
    call MPI_Comm_rank(MPI_COMM_WORLD,myrank,ierr)
    call PARCEL_Init(myrank)
    call make_matrix_DIA(n,gn,nnd,istart,diaA,offset)
    allocate(vecx(n))
    allocate(vecb(n))
    allocate(reshistory(itrmax))
    ipreflag=0
    ILU_method=1
    addL=0
    iflagAS=1
    itrmax=100
    rtolmax=1.0d-8
    iovlflag=0
    precon_thblock=-1
    independ_nvec=-1
    nblock = 2000
    vecb=1.0d0
    vecx=1.0d0
 #ifdef CUDA_SOLVER
    call cudaMalloc_FP64(cu_vecx,n)
    call cudaMalloc_FP64(cu_vecb,n)
    call cudaMalloc_FP64(cu_diaA,n*nnd)
    call cudaMalloc_PARCEL_INT(cu_offset,nnd)
    call cudaMemcpy_h2d_fp64(vecx,cu_vecx,n,offset_zero,offset_zero)
    call cudaMemcpy_h2d_fp64(vecb,cu_vecb,n,offset_zero,offset_zero)
    call cudaMemcpy_h2d_fp64(diaMat%diaA,cu_diaA,n*nnd,offset_zero,offset_zero)
    call cudaMemcpy_h2d_PARCEL_INT(diaMat%offset,cu_offset,nnd,offset_zero,offset_zero)
 #endif
    call parcel_dcg_dia(&
         MPI_COMM_WORLD, &
 #ifdef CUDA_SOLVER
         cu_vecx%mem,cu_vecb%mem, &
         n,gn,istart,&
         cu_diaA%mem,cu_offset%mem,nnd, &
 #else
         vecx,vecb,&
         n,gn,istart,&
         diaA,offset,nnd,&
 #endif
         ipreflag,ILU_method,addL,iflagAS,&
         itrmax,rtolmax,&
         reshistory,&
         iovlflag,
         precon_thblock, independ_nvec, &
         nblock, &
         iret&
         )
    call PARCEL_Free()
    call MPI_Finalize(ierr)
 #ifdef CUDA_SOLVER
    call cudaMemcpy_d2h_fp64(cu_vecx,vecx,n,offset_zero,offset_zero)
    call cudaFree(cu_vecx)
    call cudaFree(cu_vecb)
    call cudaFree(cu_diaA)
    call cudaFree(cu_offset)
 #endif
    deallocate(reshistory)
  end program main

8.3 DDM形式 (Fortran)

　DDM形式のFortranのサンプルコードを以下に示す. set_rank_xyz,set_grids_decompostion,set_neighbor_ranksはPARCEL形式の隣接プロセスを計算する任意のサブルーチン. make_matrix_DDMはPARCELのDDM形式の行列を生成する任意のサブルーチンとする.

    program main
      use mpi
      use PARCEL_integer
      use krylov
      use ddm_wrapper
      implicit none
      integer ityp_eq
      integer gnx,gny,gnz
      integer itrmax
      integer MAX_NITER
      integer nx,ny,nz
      integer nx0,ny0,nz0
      integer n,gn
      integer m
      integer m_
      integer npes,myrank
      integer ierr
      integer i
      integer ipreflag
      integer addL
      integer iflagAS
      real*8  rtolmax
      real*8  abstolmax
      integer solver,ityp_solver
      integer iret
      real*8,allocatable  :: reshistory_DDM(:)
      integer ILU_method
      integer iovlflag
      real*8,allocatable  :: vecx(:)
      real*8,allocatable  :: vecb(:)
      integer niter,restart
      integer npey,npez
      integer precon_thblock
      integer independ_nvec
      integer nblock
      integer npe_x,npe_y,npe_z
      integer div_direc_th
      integer rank_x, rank_y, rank_z
      integer,parameter :: nnd = 7
      integer,parameter :: margin = 1
      integer,parameter :: ndim=3
      integer npe_dim(ndim)
      integer gn_grids(ndim)
      integer n_grids(ndim)
      integer istart_grids(ndim)
      integer,parameter :: num_neighbor_ranks_max = 3**ndim -1
      integer neighbor_ranks( num_neighbor_ranks_max )
      integer num_neighbor_ranks
      real*8, allocatable  :: diaA(:)
      integer, allocatable :: offset(:)
      integer, allocatable :: offset_dim(:,:)
 #ifdef CUDA_SOLVER
      type(devicemem) cu_vecx
      type(devicemem) cu_vecb
      type(devicemem) cu_diaA
      type(devicemem) cu_offset
      integer(kind=intType),parameter :: offset_zero=0
 #endif
      namelist/input/ &
           ityp_eq, &
           gnx, &
           gny, &
           gnz, &
           MAX_NITER, &
           rtolmax, &
           abstolmax, &
           ipreflag, &
           independ_nvec, &
           addL, &
           iflagAS, &
           ILU_method, &
           iovlflag, &
           ityp_solver, &
           precon_thblock, &
           nblock &
      namelist/input_ddm/ &
           npe_x,&
           npe_y,&
           npe_z,&
           div_direc_th
      open(33,file='input_namelist')
      read(33,nml=input)
      read(33,nml=input_ddm)
      close(33)
      stime = 0.0d0
      etime = 0.0d0
      itrmax = MAX_NITER
      solver = ityp_solver
      call MPI_Init(ierr)
      call MPI_Comm_size(MPI_COMM_WORLD,npes,ierr)
      call MPI_Comm_rank(MPI_COMM_WORLD,myrank,ierr)
      call PARCEL_Init(myrank)
      gn_grids(1)  = gnx
      gn_grids(2)  = gny
      gn_grids(3)  = gnz
      call set_rank_xyz( myrank, npe_x, npe_y, npe_z, rank_x, rank_y, rank_z )
      call set_grids_decompostion( &
           npe_x, npe_y, npe_z, &
           rank_x, rank_y, rank_z,  &
           ndim, gn_grids, npe_dim, n_grids,  &
           istart_grids, margin, m )
      num_neighbor_ranks = -1
      call set_neighbor_ranks( &
           npe_x, npe_y, npe_z, rank_x, rank_y, rank_z, &
           num_neighbor_ranks, &
           neighbor_ranks, num_neighbor_ranks_max)
      allocate( diaA(nnd*m) )
      allocate( offset(nnd) )
      allocate( offset_dim(nnd,ndim) )
      allocate( reshistory_DDM(MAX_NITER) )
      allocate( vecx(m) )
      allocate( vecx0(m) )
      allocate( vecb(m) )
      diaA = 0.0d0
      vecx = 0.0d0
      vecx0 = 0.0d0
      vecb = 0.0d0
      reshistory_DDM = 0
      offset = 0
      offset_dim = 0
      call make_matrix_DDM( &
           gnx,gny,gnz, &
           ndim, n_grids,  &
           istart_grids, &
           margin, nnd, m, &
           offset_dim, offset, &
           diaA, vecb, vecx, vecx0,&
           rank_x, rank_y, rank_z, &
           npe_x, npe_y, npe_z )
 #ifdef CUDA_SOLVER
      call cudaMalloc_FP64(cu_vecx,m)
      call cudaMalloc_FP64(cu_vecb,m)
      call cudaMalloc_FP64(cu_diaA,m*nnd)
      call cudaMalloc_INT(cu_offset,nnd)
      call cudaMemcpy_h2d_fp64(vecx,cu_vecx,m,offset_zero,offset_zero)
      call cudaMemcpy_h2d_fp64(vecb,cu_vecb,m,offset_zero,offset_zero)
      call cudaMemcpy_h2d_fp64(diaA,cu_diaA,m*nnd,offset_zero,offset_zero)
      call cudaMemcpy_h2d_INT(offset,cu_offset, nnd,offset_zero,offset_zero)
 #endif
      call parcel_dcg_ddm(&
           MPI_COMM_WORLD, &
 #ifdef CUDA_SOLVER
           cu_vecx%mem,cu_vecb%mem, &
           m, &
           nnd,  cu_offset%mem, cu_diaA%mem, &
 #else
           vecx,vecb,  &
           m, &
           nnd, offset, diaA, &
 #endif
           num_neighbor_ranks, neighbor_ranks, margin, ndim,&
           npe_dim, n_grids, offset_dim,&
           div_direc_th, &
           ipreflag, &
           addL, iflagAS,&
           MAX_NITER, rtolmax, reshistory_DDM, &
           iovlflag, &
           precon_thblock, &
           independ_nvec, &
           nBlock, &
           iret &
           )
 #ifdef CUDA_SOLVER
      call cudaMemcpy_d2h_fp64(cu_vecx,vecx,m,offset_zero,offset_zero)
      call cudaFree(cu_vecx)
      call cudaFree(cu_vecb)
      call cudaFree(cu_diaA)
      call cudaFree(cu_offset)
 #endif
      call MPI_Barrier(MPI_COMM_WORLD,ierr)
      deallocate( diaA )
      deallocate( offset )
      deallocate( offset_dim )
      deallocate( reshistory_DDM )
      deallocate( vecx )
      deallocate( vecx0 )
      deallocate( vecb )
      call PARCEL_Free()
      call MPI_Finalize(ierr)
    end program main

9 サブルーチン使用方法(C言語)

　C言語からFortranのサブルーチンを呼ぶにはサブルーチン名の末尾に”_“を追加して引数をポインタで渡せば良い. PARCELのサブルーチンについても同様に呼び出すことで利用可能である. ただし,MPIの利用に関しては,MPIライブラリの関数MPI_Comm_c2fを用いてC用のMPIコミュニケータをFortran用のMPIコミュニケータに変換する必要がある. 以下ではC言語におけるPARCELのサブルーチン使用方法についてCG法のサンプルコードを例にして説明する. ここで,CUDA_SOLVERを有効にすればCUDA版のサンプルコードになる.

9.1 CRS形式 (C言語)

　CRS形式のサンプルコードを以下に示す. make_matrix_CRSはPARCELの行列分割に対応したCRS形式の行列を生成する任意の関数.

 
   #include < stdio.h >
   #include < stdlib.h >
   #include "mpi.h"
   #ifdef CUDA_SOLVER
   #include < cuda.h >
   #include < cuda_runtime.h >
   #endif
   int main( int argc, char *argv[]){
     int npes;
     int myrank;
     MPI_Fint icomm_fort;
     MPI_Init(&argc,&argv);
     MPI_Comm_size(MPI_COMM_WORLD, &npes);
     MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
     icomm_fort=MPI_Comm_c2f(MPI_COMM_WORLD);
     parcel_init_c_(&myrank);
     int n,gn,nnz,istart;
     double* crsA;
     int* crsRow_ptr;
     int* crsCol;
     make_matrix_CRS(n,gn,nnz,istart,crsA,crsRow_ptr,crsCol)
     double *vecx = (double *)malloc(sizeof(double)*n);
     double *vecb = (double *)malloc(sizeof(double)*n);
     for(int i=0;i < n;i++){
       vecx[i] = 0.0;
       vecb[i] = 1.0;
     }
     int ipreflag = 0;
     int ilu_method = 1;
     int addl = 0;
     int iflagas = 1;
     int max_niter = 100;
     int rtolmax = 1.0e-8;
     double *reshistory = (double *)malloc(sizeof(double)*max_niter);
     int iovlflag = 0;
     int precon_thblock = -1;
     int independ_nvec = -1;
     int nBlock = 2000;
 #ifdef CUDA_SOLVER
     typedef struct Fortran_GPUMem{
       void *a_d;
     } devicemem;
     size_t len = sizeof(double)*n;
     devicemem cu_vecx;cudaMalloc(&(cu_vecx.a_d),len);
     cudaMemcpy( cu_vecx.a_d, vecx,  len, cudaMemcpyHostToDevice);
     devicemem cu_vecb;cudaMalloc(&(cu_vecb.a_d),len);
     cudaMemcpy( cu_vecb.a_d, vecb,  len, cudaMemcpyHostToDevice);
     len = sizeof(double)*nnz;
     devicemem cu_crsa;cudaMalloc(&(cu_crsa.a_d),len);
     cudaMemcpy( cu_crsa.a_d, crsa,  len, cudaMemcpyHostToDevice);
     len = sizeof(parcel_int)*nnz;
     devicemem cu_crscol;cudaMalloc(&(cu_crscol.a_d),len);
     cudaMemcpy( cu_crscol.a_d, crscol,  len, cudaMemcpyHostToDevice);
     len = sizeof(parcel_int)*(n+1);
     devicemem cu_crsrow_ptr;cudaMalloc(&(cu_crsrow_ptr.a_d),len);
     cudaMemcpy( cu_crsrow_ptr.a_d, crsrow_ptr,  len, cudaMemcpyHostToDevice);
 #endif
     parcel_dcg_c_(
          &icomm_fort ,
 #ifdef CUDA_SOLVER
          &cu_vecx,
          &cu_vecb,
          &n,
          &gn,
          &nnz,
          &istart,
          &cu_crsa ,
          &cu_crsrow_ptr,
          &cu_crscol,
 #else
          vecx,
          vecb,
          &n,
          &gn,
          &nnz,
          &istart,
          crsa ,
          crsrow_ptr,
          crscol,
 #endif
          &ipreflag ,
          &ilu_method ,
          &addl ,
          &iflagas ,
          &max_niter ,
          &rtolmax ,
          reshistory,
          &iovlflag ,
          &precon_thblock ,
          &independ_nvec ,
          &nBlock,
          &iret
          );
 #ifdef CUDA_SOLVER
     len = sizeof(double)*n;
     cudaMemcpy(  vecx,cu_vecx.a_d,  len, cudaMemcpyDeviceToHost);
     cudaFree(cu_vecx.a_d);
     cudaFree(cu_vecb.a_d);
     cudaFree(cu_crsa.a_d);
     cudaFree(cu_crscol.a_d);
     cudaFree(cu_crsrow_ptr.a_d);
 #endif
     free(vecx);
     free(vecb);
     free(reshistory);
     parcel_free_c_();
     MPI_Finalize();
     }

9.2 DIA形式 (C言語)

　DIA形式のサンプルコードを以下に示す. make_matrix_DIAはPARCELのブロック分割に対応したDIA形式の行列を生成する任意の関数.

 
   #include < stdio.h >
   #include < stdlib.h >
   #include "mpi.h"
   #ifdef CUDA_SOLVER
   #include < cuda.h >
   #include < cuda_runtime.h >
   #endif
   int main( int argc, char *argv[]){
     int npes;
     int myrank;
     MPI_Fint icomm_fort;
     MPI_Init(&argc,&argv);
     MPI_Comm_size(MPI_COMM_WORLD, &npes);
     MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
     icomm_fort=MPI_Comm_c2f(MPI_COMM_WORLD);
     parcel_init_c_(&myrank);
     int n,gn,nnd,istart;
     double* diaA;
     int* offset;
     make_matrix_DIA(n,gn,nnd,istart,diaA,offset)
     double *vecx = (double *)malloc(sizeof(double)*n);
     double *vecb = (double *)malloc(sizeof(double)*n);
     for(int i=0;i < n;i++){
       vecx[i] = 0.0;
       vecb[i] = 1.0;
     }
     int ipreflag = 0;
     int ilu_method = 1;
     int addl = 0;
     int iflagas = 1;
     int max_niter = 100;
     int rtolmax = 1.0e-8;
     double *reshistory = (double *)malloc(sizeof(double)*max_niter);
     int iovlflag = 0;
     int precon_thblock = -1;
     int independ_nvec = -1;
     int nBlock = 2000;
 #ifdef CUDA_SOLVER
     size_t len = sizeof(double)*n;
     cudaMemcpy( cu_vecx.a_d, vecx,  len, cudaMemcpyHostToDevice);
     cudaMemcpy( cu_vecb.a_d, vecb,  len, cudaMemcpyHostToDevice);
     len = sizeof(double)*n*nnd;
     devicemem cu_diaA;cudaMalloc(&(cu_diaA.a_d),len);
     cudaMemcpy( cu_diaA.a_d, diaA,  len, cudaMemcpyHostToDevice);
     len = sizeof(parcel_int)*nnd;
     devicemem cu_offset;cudaMalloc(&(cu_offset.a_d),len);
     cudaMemcpy( cu_offset.a_d, offset,  len, cudaMemcpyHostToDevice);
 #endif
     parcel_dcg_dia_c_(
          &icomm_fort ,
 #ifdef CUDA_SOLVER
          &cu_vecx,
          &cu_vecb,
          &n,
          &gn,
          &istart,
          &cu_diaA,
          &cu_offset,
 #else
          vecx,
          vecb,
          &n,
          &gn,
          &istart,
          diaA ,
          offset,
 #endif
          &nnd,
          &ipreflag ,
          &ilu_method ,
          &addl ,
          &iflagas ,
          &max_niter ,
          &rtolmax ,
          reshistory,
          &iovlflag ,
          &precon_thblock ,
          &independ_nvec ,
          &nBlock,
          &iret
          );
 #ifdef CUDA_SOLVER
     len = sizeof(double)*n;
     cudaMemcpy(  vecx,cu_vecx.a_d,  len, cudaMemcpyDeviceToHost);
     cudaFree(cu_vecx.a_d);
     cudaFree(cu_vecb.a_d);
     cudaFree(cu_diaA.a_d);
     cudaFree(cu_offset.a_d);
 #endif
     free(vecx);
     free(vecb);
     free(reshistory);
     parcel_free_c_();
     MPI_Finalize();
     }

9.3 DDM形式 (C言語)

　DDM形式のC言語のサンプルコードを以下に示す. set_rank_xyz,set_grids_decompostion,set_neighbor_ranksはPARCEL形式の隣接プロセスを計算する任意の関数. make_matrix_ddmはPARCELのDDM形式の行列を生成する任意の関数とする.

 
   #include < stdio.h >
   #include < stdlib.h >
   #include < math.h >
   #include "mpi.h"
   #ifdef CUDA_SOLVER
   #include < cuda.h >
   #include < cuda_runtime.h >
   #endif
   int main( int argc, char *argv[]){
     int ityp_eq;
     int itrmax;
     int max_niter;
     int nx,ny,nz;
     int nx0,ny0,nz0;
     int n,gn;
     int istart,iend;
     int npes,myrank;
     int ierr;
     int i;
     int iret;
     int myranky,myrankz;
     int npey,npez;
     MPI_Fint icomm_fort;
     int gnx = 100;
     int gny = 100;
     int gnz = 100;
     int max_niter = 100;
     int rtolmax = 1.0e-8;
     int ipreflag = 0;
     int independ_nvec = -1;
     int addl = 0;
     int iflagas = 1;
     int ilu_method = 1;
     int iovlflag = 0;
     int precon_thblock = -1;
     int nBlock = 2000;
     int npe_x = 1;
     int npe_y = 4;
     int npe_z = 1;
     int div_direc_th = 1;
     MPI_Init(&argc,&argv);
     MPI_Comm_size(MPI_COMM_WORLD, &npes);
     MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
     icomm_fort=MPI_Comm_c2f(MPI_COMM_WORLD);
     parcel_init_c_(&myrank);
     int nnd = 7;
     int rank_x,rank_y,rank_z;
     set_rank_xyz(
          &myrank,
          &npe_x,
          &npe_y,
          &npe_z,
          &rank_x,
          &rank_y,
          &rank_z);
     int ndim=3 ;
     int margin = 1;
     int *npe_dim = (int *)malloc(sizeof(int)*ndim);
     int *gn_grids = (int *)malloc(sizeof(int)*ndim);
     int *n_grids = (int *)malloc(sizeof(int)*ndim);
     int *istart_grids = (int *)malloc(sizeof(int)*ndim);
     for(i=0;i < ndim;i++){
       npe_dim[i] = 0;
       gn_grids[i] = 0;
       n_grids[i] = 0;
       istart_grids[i] = 0;
     }
     int m=-1;
     gn_grids[0]  = gnx;
     gn_grids[1]  = gny;
     gn_grids[2]  = gnz;
     set_grids_decompostion(
                &npe_x, &npe_y, &npe_z,
                &rank_x, &rank_y, &rank_z,
                &ndim, gn_grids, npe_dim, n_grids,
                istart_grids, &margin, &m );
     int num_neighbor_ranks_max = pow(3,ndim) -1;
     int *neighbor_ranks = (int *)malloc(sizeof(int)*num_neighbor_ranks_max);
     int num_neighbor_ranks=-1;
     for(i=0;i < num_neighbor_ranks_max;i++){
       neighbor_ranks[i] = -1;
     }
     set_neighbor_ranks(
            &npe_x,
            &npe_y,
            &npe_z,
            &rank_x,
            &rank_y,
            &rank_z,
            &num_neighbor_ranks,
            neighbor_ranks,
            &num_neighbor_ranks_max);
     int *offset = (int *)malloc(sizeof(int)*nnd);
     for(i=0;i < nnd;i++){
       offset[i] = 0;
     }
     int *offset_dim = (int *)malloc(sizeof(int)*nnd*ndim);
     for(i=0;i < nnd*ndim;i++){
       offset_dim[i] = 0;
     }
     double *diaA  = (double *)malloc(sizeof(double)*m*nnd);
     for(i=0;i < m*nnd;i++){
       diaA[i] = 0.0;
     }
     double *vecx  = (double *)malloc(sizeof(double)*m);
     double *vecb  = (double *)malloc(sizeof(double)*m);
     for(i=0;i < m;i++){
       vecx[i] = 0.0;
       vecb[i] = 0.0;
     }
     double *reshistory_ddm  = (double *)malloc(sizeof(double)*max_niter);
     for(i=0;i < max_niter;i++){
       reshistory_ddm[i] = 0.0;
     }
     int nnd_ = nnd;
     make_matrix_ddm(
             &gnx,
             &gny,
             &gnz,
             &ndim,
             n_grids,
             istart_grids,
             &margin,
             &nnd_,
             &m,
             offset_dim,
             offset,
             diaA,
             vecb,
             vecx,
             vecx0,
             &rank_x,
             &rank_y,
             &rank_z,
             &npe_x,
             &npe_y,
             &npe_z);
 #ifdef CUDA_SOLVER
     typedef struct Fortran_GPUMem{
       void *a_d;
     } devicemem;
     size_t len = sizeof(double)*m;
     devicemem cu_vecx;cudaMalloc(&(cu_vecx.a_d),len);
     cudaMemcpy( cu_vecx.a_d, vecx,  len, cudaMemcpyHostToDevice);
     devicemem cu_vecb;cudaMalloc(&(cu_vecb.a_d),len);
     cudaMemcpy( cu_vecb.a_d, vecb,  len, cudaMemcpyHostToDevice);
     len = sizeof(double)*m*nnd;
     devicemem cu_diaA;cudaMalloc(&(cu_diaA.a_d),len);
     cudaMemcpy( cu_diaA.a_d, diaA,  len, cudaMemcpyHostToDevice);
     len = sizeof(int)*nnd;
     devicemem cu_offset;cudaMalloc(&(cu_offset.a_d),len);
     cudaMemcpy( cu_offset.a_d, offset,  len, cudaMemcpyHostToDevice);
 #endif
     parcel_dcg_ddm_c_(
              &icomm_fort,
 #ifdef CUDA_SOLVER
              &cu_vecx,
              &cu_vecb,
              &m,
              &nnd,
              &cu_offset,
              &cu_diaA,
 #else
              vecx,
              vecb,
              &m,
              &nnd,
              offset,
              diaA,
 #endif
              &num_neighbor_ranks,
              neighbor_ranks,
              &margin,
              &ndim,
              npe_dim,
              n_grids,
              offset_dim,
              &div_direc_th,
              &ipreflag,
              &addl,
              &iflagas,
              &max_niter,
              &rtolmax,
              reshistory_ddm,
              &iovlflag,
              &precon_thblock,
              &independ_nvec ,
              &nBlock ,
              &iret
              );
 #ifdef CUDA_SOLVER
     len = sizeof(double)*m;
     cudaMemcpy(  vecx,cu_vecx.a_d,  len, cudaMemcpyDeviceToHost);
     cudaFree(cu_diaA.a_d);
     cudaFree(cu_offset.a_d);
     cudaFree(cu_vecx.a_d);
     cudaFree(cu_vecb.a_d);
 #endif
     free(offset);
     free(offset_dim);
     free(diaA);
     free(vecx);
     free(vecb);
     free(reshistory_ddm);
     free(neighbor_ranks);
     free(npe_dim);
     free(gn_grids);
     free(n_grids);
     free(istart_grids);
     MPI_Barrier(MPI_COMM_WORLD);
     parcel_free_c_();
     MPI_Finalize();
     return 0;
   }

10 参考文献

[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] Y. Saad, “Iterative methods for sparse linear systems”, SIAM (2003)
[3] M. Hoemmen, “Communication-avoiding Krylov subspace methods”. Ph.D.thesis, University of California, Berkeley（2010）
[4] 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).
[5] R. Suda, L. Cong, D. Watanabe, et al., “Communication-Avoiding CG Method: New Direction of Krylov Subspace Methods towards Exa-scale Computing”, RIMS Kokyuroku ,pp. 102-111, (2016).
[6] Y. Idomura, T. Ina, S. Yamashita, et al., “Communication avoiding multigrid preconditioned conjugate gradient method for extreme scale multiphase CFD simulations”. ScalA 18:9th Workshop on Latest Advances in Scalable Algorithms for Large Scale Systems,pp. 17-24.(2018)
[7] T. Ina, Y. Idomura, T. Imamura, et al., “Iterative methods with mixed-precision preconditioning for ill-conditioned linear systems in multiphase CFD simulations”, ScalA21:12th Workshop on Latest Advances in Scalable Algorithms for Large Scale Systems .(2021)
[8] A. Stathopoulos, K. Wu, “A block orthogonalization procedure with constant synchronization requirements”. SIAM J. Sci. Comput. 23, 2165–2182 (2002)
[9] T. Fukaya, Y. Nakatsukasa, Y. Yanagisawa, et al., “CholeskyQR2: A Simple and Communication-Avoiding Algorithm for Computing a Tall-Skinny QR Factorization on a Large-Scale Parallel System,” ScalA 14:5th Workshop on Latest Advances in Scalable Algorithms for Large Scale Systems, pp. 31-38,(2014)