1 クリロフ部分空間法

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

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 前処理付安定化双共役勾配法

　安定化双共役勾配法(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 前処理付一般化最小残差法

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

　通信削減型一般化最小残差法(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 前処理付チェビシェフ基底共役勾配法法

　チェビシェフ基底共役勾配法(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{)}$

2 前処理

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

2.1 ポイントヤコビ前処理

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

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

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

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

　係数行列$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 ブロックヤコビ前処理

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

2.5 加法的シュワルツ前処理

　前処理行列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$を考慮しない結果を採用する.

3 データ構造

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

3.1 Compressed Row 形式

　Compressed Row 形式(圧縮行形式)は,非ゼロ要素を行方向に圧縮した後,列方向にデータを並べて格納する方法である. 非ゼロ要素の列位置・値ともに全体を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\}$

3.2 Diagonal形式

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

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

4.1 ディレクトリ構成

　展開したアーカイブの内容を図6に示す.

4.2 コンパイル方法

　PARCELのコンパイルにはCコンパイラ及びプリプロセッサ対応のFORTRANコンパイラ. MPI,LAPACK,OpenMPライブラリを必要とする. gfortranでコンパイルする場合はコンパイルオプションに-cppを追加することでプリプロセッサが適用される. コンパイルオプションはarch/make_configを編集することで変更する. make_configの記述例を図7に示す. コンパイル方法は展開したアーカイブのトップディレクトリでmakeコマンドを実行する. 再コンパイルする場合はmake cleanを実行後にmakeコマンドを実行する. 各ディレクトリのexampleファイルに実行ファイルが生成されていればコンパイル成功である. 20億次元以上の大規模問題では64ビット整数が必要となる. 例えば,本パッケージに格納されている大規模テストexample5のコンパイルにはFDEFINEの指定が必要である。

4.3 使用方法

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

5 PARCELサブルーチン（直接インターフェース）

　PARCELに整備されているサブルーチンの直接呼び出しインターフェースについて説明する.

5.1 parcel_dcg

• 呼び出し方法

call parcel_dcg( icomm, vecx,vecb, n,gn,nnz,istart, crsA,crsRow_ptr,crsCol, ipreflag,ILU_method,addL,iflagAS, itrmax,rtolmax, reshistory, iovlflag,iret )

• 機能

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

5.2 parcel_dcg_dia

• 呼び出し方法

call parcel_dcg_dia( icomm, vecx,vecb, n,gn,istart, diaA,offset,nnd, ipreflag,ILU_method,addL,iflagAS, itrmax,rtolmax, reshistory, iovlflag,iret )

• 機能

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

5.3 parcel_ddcg

• 呼び出し方法

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,iret, PRECISION_A, PRECISION_b, PRECISION_x, PRECISION_PRECONDITION )

• 機能

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

5.4 parcel_ddcg_dia

• 呼び出し方法

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,iret, PRECISION_A, PRECISION_b, PRECISION_x, PRECISION_PRECONDITION )

• 機能

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

5.5 parcel_dbicgstab

• 呼び出し方法

call parcel_dbicgstab( icomm, vecx,vecb, n,gn,nnz,istart, crsA,crsRow_ptr,crsCol, ipreflag,ILU_method,addL,iflagAS, itrmax,rtolmax, reshistory, iovlflag,iret )

• 機能

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

*加法シュワルツを使用する際はILU分解の手法をD-ILUから選択する必要がある。＊

5.6 parcel_dbicgstab_dia

• 呼び出し方法

call parcel_dbicgstab_dia( icomm, vecx,vecb, n,gn,istart, diaA,offset,nnd, ipreflag,ILU_method,addL,iflagAS, itrmax,rtolmax, reshistory, iovlflag,iret )

• 機能

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

5.7 parcel_ddbicgstab

• 呼び出し方法

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,iret, PRECISION_A, PRECISION_b, PRECISION_x, PRECISION_PRECONDITION )

• 機能

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

5.8 parcel_ddbicgstab_dia

• 呼び出し方法

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,iret, PRECISION_A, PRECISION_b, PRECISION_x, PRECISION_PRECONDITION )

• 機能

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