1 クリロフ部分空間法

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


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

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

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


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

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

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


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

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

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

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


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

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

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

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

  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\ \ }}\)

\(\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}\)

 \(\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}\)

  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 前処理付チェビシェフ基底共役勾配法[4,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}}\)

\(\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 通信削減型アーノルディ法[2]

 通信削減型アーノルディ法(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]

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


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

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

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

 係数行列\(A\)をその対角行列\(D\),下三角行列\(L\),上三角行列\(U\)を用いて \[A = L_{A} + D_{A} + U_{A}\] と分解する.この\(L_{A}\),\(U_{A}\)および対角行列\(D\) (\(D \neq D_{A}\))を用いて,前処理行列\(K\)\[K = \left( D + L_{A} \right)D^{- 1}\left( D + U_{A} \right)\] の形に決定する.ここでは,対角行列\(D\)の決定方法として,次の2とおりの条件のどちらを選択する.

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

となる.


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

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


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

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

図1 加法的シュワルツ前処理

BASIC

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

図2 BASIC

RESTRICT

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

図3 RESTRICT

INTERPOLATE

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

図4 INTERPOLATE

NONE

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

図5 NONE


2.6 細分化前処理[6]

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

 ブロックヤコビ前処理 

BlockJacobi

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

Subdividing BlockJacobi
                図6細 分化前処理


3 QR分解

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


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

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

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


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

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

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


3.3 Tall Skinny QR (TSQR)法[2]

 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
SequentialTSQR       ParallelTSQR


3.4 コレスキーQR法[7]

 コレスキー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法[8]

 コレスキー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プロセスの場合

行列の非ゼロ要素  :: \(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プロセスの場合

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

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

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

行列の非ゼロ要素  :: \(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プロセスの場合

対角要素    :: \(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プロセスの場合

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

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

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

対角要素    :: \(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次元の構造格子で与えられる行列の並列処理に必要となる通信を最小化することができる. 関係するパラメータ及び指定法を以下に説明する.

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

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


ixmax = 0; iymax = 0; izmax = 0; num = 0;

neighbor_ranks = -1

if ( npe_x > 1 ) ixmax=1

if ( npe_y > 1 ) iymax=1

if ( npe_z > 1 ) izmax=1

do iz=-izmax, izmax

  do iy=-iymax, iymax

    do ix=-ixmax, ixmax

      n1x = rank_x +ix; n1y = rank_y +iy; n1z = rank_z +iz;

      if ( ix*ix +iy*iy + iz*iz == 0 ) cycle

      num = num +1

      if ( n1z >=0 .and. n1z < npe_z ) then

        if ( n1y >=0 .and. n1y < npe_y ) then

          if ( n1x >=0 .and. n1x < npe_x ) then

            iirank = npe_x *npe_y *n1z +npe_x *n1y +n1x;

            neighbor_ranks(num) = iirank;

          endif

        endif

      endif

    end do

  end do

end do


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

i  ioff_grids(*,1)  ioff_grids(*,2)  ioff_grids(*,3)
 基準 1 0 0 0
 東  2 1 0 0
 西  3 -1 0 0
 南  4 0 1 0
 北  5 0 -1 0
 上  6 0 0 1
 下  7 0 0 -1

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プロセスの場合

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

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

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

ステンシル成分 ::

i  ioff_grids(*,1)
 基準 1 0
 東  2 1
 西  3 -1

2プロセスの場合

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

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

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

ステンシル成分 ::

i  ioff_grids(*,1)
 基準 1 0
 東  2 1
 西  3 -1

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

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

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

ステンシル成分 ::

i  ioff_grids(*,1)
 基準 1 0
 東  2 1
 西  3 -1


5 インストール

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

5.1 コンパイル方法

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

  CC     = mpiicc
  FC     = mpiifort
 
  CFLAGS  = -O3
  FFLAGS  = -O3 -fpp -qopenmp -xHost  -mcmodel=large
 
  FP_MODE = -fp-model precise
 
  INCLUDEDIR =
  LD_MPI     = -lmpifort -lifcore
  LD_LAPACK  = -mkl
 
  #FDEFINE = -DPARCEL_INT8
  INSTALLDIR = ../PARCEL_1.1

    オプション名     説明
CC Cコンパイラコマンド名
FC FORTRANコンパイラコマンド名
CFLAGS Cコンパイラコンパイルオプション
FFLAGS FORTRANコンパイラコンパイルオプション
FP_MODE 演算順序変更を伴う最適化を抑制するオプション(4倍精度ライブラリのみに適用される)
INCLUDEDIR インクルードディレクトリ
LD_MPI MPIライブラリ
LD_LAPACK LAPACKライブラリ
FDEFINE -DPARCEL_INT8を指定することで64ビット整数型に対応.
INSTALLDIR PARCELのインストール先ディレクトリ


5.2 使用方法

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


6 PARCELサブルーチン

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

 parcel_dから始まる関数名は倍精度、parcel_ddから始まる関数名は4倍精度(Double-Double形式)で実装されている.

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

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

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

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

 i番目の固有値が実数であれば対応する固有ベクトルはEvec(n*(i-1)+1:n*i)に格納される.

 i番目の固有値が複素数であり,i+1番目の固有値がi番目の固有値の共役複素数の場合のi,i+1番目の固有ベクトルは次のように格納される.

固有値 実部 虚部
i番目     Evec(n*(i-1)+1:n*i)         Evec(n*i+1:n*(i+1))    
(i+1)番目     Evec(n*(i-1)+1:n*i)        -Evec(n*i+1:n*(i+1))    


6.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, precon_thblock, independ_nvec, nBlock, iret )


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

引数名(次元) 入力/出力 説明
icomm integer in MPIコミュニケータ
vecx(n) double precision in/out in:反復の初期値,out:解ベクトル
vecb(n) double precision in 右辺ベクトル
n integer*4 / integer*8 in 各プロセスが担当するベクトルの大きさ
gn integer*4 / integer*8 in 全体のベクトルの大きさ
nnz integer*4 / integer*8 in 各プロセスが担当する係数行列の非零要素数
istart integer*4 / integer*8 in 各プロセスが担当する行列の開始行
crsA(nnz) double precision in CRS 形式に格納された行列の非零要素
crsRow_ptr(n+1) integer*4 / integer*8 in CRS 形式のポインタテーブル
crsCol(nnz) integer*4 / integer*8 in CRS 形式に格納された分割行列の非零要素の列番号
ipreflag integer in 前処理指定フラグ(0:無し,1:ポイントヤコビ,2:ブロックヤコビ,3:加法シュワルツ)
ilu_method integer in 不完全LU分解指定フラグ(0:ILU(0),1:D-ILU(対角成分が一致),2:D-ILU(行内の要素和が一致),加法シュワルツは1,2のみ利用可能)
addL integer in 加法シュワルツの重なり幅
iflagAS integer in 加法シュワルツ指定フラグ(1:BASIC,2:RESTRICT,3:INTERPOLATE,4:NONE)
itrmax integer in/out in:最大反復回数,out:反復回数
rtolmax double precision in 収束判定条件(相対残差ノルム)
reshistory(itrmax) double precision out 残差履歴
iovlflag integer in 通信隠蔽フラグ(0:無し,1:全プロセス,2:各プロセス)
precon_thblock integer in ブロックヤコビ前処理,加法シュワルツ前処理のブロック数
independ_nvec integer in 細分化前処理で独立に計算する対角ブロックのサイズ
nBlock integer in サイクリックにスレッドに割り当てる行数
iret integer out エラーコード(0:正常終了)


6.2 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, precon_thblock, independ_nvec, nBlock, iret )


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

引数名(次元) 入力/出力 説明
icomm integer in MPIコミュニケータ
vecx(n) double precision in/out in:反復の初期値,out:解ベクトル
vecb(n) double precision in 右辺ベクトル
n integer*4 / integer*8 in 各プロセスが担当するベクトルの大きさ
gn integer*4 / integer*8 in 全体のベクトルの大きさ
nnz integer*4 / integer*8 in 各プロセスが担当する係数行列の非零要素数
istart integer*4 / integer*8 in 各プロセスが担当する行列の開始行
crsA(nnz) double precision in CRS 形式に格納された行列の非零要素
crsRow_ptr(n+1) integer*4 / integer*8 in CRS 形式のポインタテーブル
crsCol(nnz) integer*4 / integer*8 in CRS 形式に格納された分割行列の非零要素の列番号
ipreflag integer in 前処理指定フラグ(0:無し,1:ポイントヤコビ,2:ブロックヤコビ,3:加法シュワルツ)
ilu_method integer in 不完全LU分解指定フラグ(0:ILU(0),1:D-ILU(対角成分が一致),2:D-ILU(行内の要素和が一致),加法シュワルツは1,2のみ利用可能)
addL integer in 加法シュワルツの重なり幅
iflagAS integer in 加法シュワルツ指定フラグ(1:BASIC,2:RESTRICT,3:INTERPOLATE,4:NONE)
itrmax integer in/out in:最大反復回数,out:反復回数
rtolmax double precision in 収束判定条件(相対残差ノルム)
reshistory(itrmax) double precision out 残差履歴
iovlflag integer in 通信隠蔽フラグ(0:無し,1:全プロセス,2:各プロセス)
precon_thblock integer in ブロックヤコビ前処理,加法シュワルツ前処理のブロック数
independ_nvec integer in 細分化前処理で独立に計算する対角ブロックのサイズ
nBlock integer in サイクリックにスレッドに割り当てる行数
iret integer out エラーコード(0:正常終了)


6.3 parcel_dgmres


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, independ_nvec, nBlock, gmres_m,gmres_GSflag, iret )


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

引数名(次元) 入力/出力 説明
icomm integer in MPIコミュニケータ
vecx(n) double precision in/out in:反復の初期値,out:解ベクトル
vecb(n) double precision in 右辺ベクトル
n integer*4 / integer*8 in 各プロセスが担当するベクトルの大きさ
gn integer*4 / integer*8 in 全体のベクトルの大きさ
nnz integer*4 / integer*8 in 各プロセスが担当する係数行列の非零要素数
istart integer*4 / integer*8 in 各プロセスが担当する行列の開始行
crsA(nnz) double precision in CRS 形式に格納された行列の非零要素
crsRow_ptr(n+1) integer*4 / integer*8 in CRS 形式のポインタテーブル
crsCol(nnz) integer*4 / integer*8 in CRS 形式に格納された分割行列の非零要素の列番号
ipreflag integer in 前処理指定フラグ(0:無し,1:ポイントヤコビ,2:ブロックヤコビ,3:加法シュワルツ)
ilu_method integer in 不完全LU分解指定フラグ(0:ILU(0),1:D-ILU(対角成分が一致),2:D-ILU(行内の要素和が一致),加法シュワルツは1,2のみ利用可能)
addL integer in 加法シュワルツの重なり幅
iflagAS integer in 加法シュワルツ指定フラグ(1:BASIC,2:RESTRICT,3:INTERPOLATE,4:NONE)
itrmax integer in/out in:最大反復回数,out:反復回数
rtolmax double precision in 収束判定条件(相対残差ノルム)
reshistory(itrmax) double precision out 残差履歴
iovlflag integer in 通信隠蔽フラグ(0:無し,1:全プロセス,2:各プロセス)
precon_thblock integer in ブロックヤコビ前処理,加法シュワルツ前処理のブロック数
independ_nvec integer in 細分化前処理で独立に計算する対角ブロックのサイズ
nBlock integer in サイクリックにスレッドに割り当てる行数
gmres_m integer in リスタートするまでの反復回数
gmres_GSflag integer in 直交化フラグ(1:MGS,2:CGS)
iret integer out エラーコード(0:正常終了)


6.4 parcel_dcagmres


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, independ_nvec, nBlock, cagmres_sstep, cagmres_tstep, cagmres_basis, cagmres_QRflag, iret )


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

引数名(次元) 入力/出力 説明
icomm integer in MPIコミュニケータ
vecx(n) double precision in/out in:反復の初期値,out:解ベクトル
vecb(n) double precision in 右辺ベクトル
n integer*4 / integer*8 in 各プロセスが担当するベクトルの大きさ
gn integer*4 / integer*8 in 全体のベクトルの大きさ
nnz integer*4 / integer*8 in 各プロセスが担当する係数行列の非零要素数
istart integer*4 / integer*8 in 各プロセスが担当する行列の開始行
crsA(nnz) double precision in CRS 形式に格納された行列の非零要素
crsRow_ptr(n+1) integer*4 / integer*8 in CRS 形式のポインタテーブル
crsCol(nnz) integer*4 / integer*8 in CRS 形式に格納された分割行列の非零要素の列番号
ipreflag integer in 前処理指定フラグ(0:無し,1:ポイントヤコビ,2:ブロックヤコビ,3:加法シュワルツ)
ilu_method integer in 不完全LU分解指定フラグ(0:ILU(0),1:D-ILU(対角成分が一致),2:D-ILU(行内の要素和が一致),加法シュワルツは1,2のみ利用可能)
addL integer in 加法シュワルツの重なり幅
iflagAS integer in 加法シュワルツ指定フラグ(1:BASIC,2:RESTRICT,3:INTERPOLATE,4:NONE)
itrmax integer in/out in:最大反復回数,out:反復回数
rtolmax double precision in 収束判定条件(相対残差ノルム)
reshistory(itrmax) double precision out 残差履歴
iovlflag integer in 通信隠蔽フラグ(0:無し,1:全プロセス,2:各プロセス)
precon_thblock integer in ブロックヤコビ前処理,加法シュワルツ前処理のブロック数
independ_nvec integer in 細分化前処理で独立に計算する対角ブロックのサイズ
nBlock integer in サイクリックにスレッドに割り当てる行数
cagmres_sstep integer in CA-GMRES法の省通信ステップ数
cagmres_tstep integer in CA-GMRES法のリスタートステップ数
cagmres_basis integer in CA-GMRES法の基底生成フラグ(0:単基底,1:ニュートン基底)
cagmres_QRflag integer in CA-GMRES法のQR分解フラグ(1:MGS,2:CGS,3:TSQR,4:CholeskyQR,5:CholeskyQR2)
iret integer out エラーコード(0:正常終了)


6.5 parcel_dcbcg


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, independ_nvec, nBlock, cbcg_kstep,cbcg_Eigenflag,power_method_itrmax, caarnoldi_sstep,caarnoldi_tstep, caarnoldi_basis,caarnoldi_QRflag, iret )


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

引数名(次元) 入力/出力 説明
icomm integer in MPIコミュニケータ
vecx(n) double precision in/out in:反復の初期値,out:解ベクトル
vecb(n) double precision in 右辺ベクトル
n integer*4 / integer*8 in 各プロセスが担当するベクトルの大きさ
gn integer*4 / integer*8 in 全体のベクトルの大きさ
nnz integer*4 / integer*8 in 各プロセスが担当する係数行列の非零要素数
istart integer*4 / integer*8 in 各プロセスが担当する行列の開始行
crsA(nnz) double precision in CRS 形式に格納された行列の非零要素
crsRow_ptr(n+1) integer*4 / integer*8 in CRS 形式のポインタテーブル
crsCol(nnz) integer*4 / integer*8 in CRS 形式に格納された分割行列の非零要素の列番号
ipreflag integer in 前処理指定フラグ(0:無し,1:ポイントヤコビ,2:ブロックヤコビ,3:加法シュワルツ)
ilu_method integer in 不完全LU分解指定フラグ(0:ILU(0),1:D-ILU(対角成分が一致),2:D-ILU(行内の要素和が一致),加法シュワルツは1,2のみ利用可能)
addL integer in 加法シュワルツの重なり幅
iflagAS integer in 加法シュワルツ指定フラグ(1:BASIC,2:RESTRICT,3:INTERPOLATE,4:NONE)
itrmax integer in/out in:最大反復回数,out:反復回数
rtolmax double precision in 収束判定条件(相対残差ノルム)
reshistory(itrmax) double precision out 残差履歴
iovlflag integer in 通信隠蔽フラグ(0:無し,1:全プロセス,2:各プロセス)
precon_thblock integer in ブロックヤコビ前処理,加法シュワルツ前処理のブロック数
independ_nvec integer in 細分化前処理で独立に計算する対角ブロックのサイズ
nBlock integer in サイクリックにスレッドに割り当てる行数
cbcg_kstep integer in CBCG法の省通信ステップ
cbcg_Eigenflag integer in 固有値計算フラグ(1:べき乗法,2:CA-ARNOLDI)
power_method_itrmax integer in べき乗法の反復回数
caarnoldi_sstep integer in CA-ARNOLDI法の省通信ステップ数
caarnoldi_tstep integer in CA-ARNOLDI法のリスタートステップ数
caarnoldi_basis integer in CA-ARNOLDI法の基底生成フラグ(0:単基底,1:ニュートン基底)
caarnoldi_QRflag integer in CA-ARNOLDI法のQR分解フラグ(1:MGS,2:CGS,3:TSQR,4:CholeskyQR,5:CholeskyQR2)
iret integer out エラーコード(0:正常終了)


6.6 parcel_dcaarnoldi


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形式で格納.

引数名(次元) 入力/出力 説明
icomm integer in MPIコミュニケータ
vecx(n) double precision in/out in:反復の初期値,out:解ベクトル
vecb(n) double precision in 右辺ベクトル
n integer*4 / integer*8 in 各プロセスが担当するベクトルの大きさ
gn integer*4 / integer*8 in 全体のベクトルの大きさ
nnz integer*4 / integer*8 in 各プロセスが担当する係数行列の非零要素数
istart integer*4 / integer*8 in 各プロセスが担当する行列の開始行
crsA(nnz) double precision in CRS 形式に格納された行列の非零要素
crsRow_ptr(n+1) integer*4 / integer*8 in CRS 形式のポインタテーブル
crsCol(nnz) integer*4 / integer*8 in CRS 形式に格納された分割行列の非零要素の列番号
ipreflag integer in 前処理指定フラグ(0:無し,1:ポイントヤコビ,2:ブロックヤコビ,3:加法シュワルツ)
ilu_method integer in 不完全LU分解指定フラグ(0:ILU(0),1:D-ILU(対角成分が一致),2:D-ILU(行内の要素和が一致),加法シュワルツは1,2のみ利用可能)
addL integer in 加法シュワルツの重なり幅
iflagAS integer in 加法シュワルツ指定フラグ(1:BASIC,2:RESTRICT,3:INTERPOLATE,4:NONE)
itrmax integer in/out in:最大反復回数,out:反復回数
iovlflag integer in 通信隠蔽フラグ(0:無し,1:全プロセス,2:各プロセス)
precon_thblock integer in ブロックヤコビ前処理,加法シュワルツ前処理のブロック数
independ_nvec integer in 細分化前処理で独立に計算する対角ブロックのサイズ
nBlock integer in サイクリックにスレッドに割り当てる行数
caarnoldi_sstep integer in CA-ARNOLDI法の省通信ステップ数
caarnoldi_tstep integer in CA-ARNOLDI法のリスタートステップ数
caarnoldi_basis integer in CA-ARNOLDI法の基底生成フラグ(0:単基底,1:ニュートン基底)
caarnoldi_QRflag integer in CA-ARNOLDI法のQR分解フラグ(1:MGS,2:CGS,3:TSQR,4:CholeskyQR,5:CholeskyQR2)
Evalr(caarnoldi_sstep*caarnoldi_tstep) double precision out 固有値の実数部
Evali(caarnoldi_sstep*caarnoldi_tstep) double precision out 固有値の虚数部
Evec(n*caarnoldi_sstep*caarnoldi_tstep) double precision out 固有ベクトル
Eerr(caarnoldi_sstep*caarnoldi_tstep) double precision out 固有値及び固有ベクトルのエラーノルム \({\left\| Ax - \lambda x \right\|}/{\left\| \text{λx} \right\|}\)
EmaxID integer out 最大固有値のID番号
EminID integer out 最小固有値のID番号
iret integer out エラーコード(0:正常終了)


6.7 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, precon_thblock, independ_nvec, nBlock, iret )


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

引数名(次元) 入力/出力 説明
icomm integer in MPIコミュニケータ
vecx(n) double precision in/out in:反復の初期値,out:解ベクトル
vecb(n) double precision in 右辺ベクトル
n integer*4 / integer*8 in 各プロセスが担当するベクトルの大きさ
gn integer*4 / integer*8 in 全体のベクトルの大きさ
istart integer*4 / integer*8 in 各プロセスが担当する行列の開始行
diaA(n*nnd) double precision in DIA形式に格納された行列の非零要素
offset(nnd) integer*4 / integer*8 in DIA形式に格納された行列の対角からのオフセット
nnd integer*4 / integer*8 in DIA形式に格納された行列の対角要素の本数
ipreflag integer in 前処理指定フラグ(0:無し,1:ポイントヤコビ,2:ブロックヤコビ,3:加法シュワルツ)
ilu_method integer in 不完全LU分解指定フラグ(0:ILU(0),1:D-ILU(対角成分が一致),2:D-ILU(行内の要素和が一致),加法シュワルツは1,2のみ利用可能)
addL integer in 加法シュワルツの重なり幅
iflagAS integer in 加法シュワルツ指定フラグ(1:BASIC,2:RESTRICT,3:INTERPOLATE,4:NONE)
itrmax integer in/out in:最大反復回数,out:反復回数
rtolmax double precision in 収束判定条件(相対残差ノルム)
reshistory(itrmax) double precision out 残差履歴
iovlflag integer in 通信隠蔽フラグ(0:無し,1:全プロセス,2:各プロセス)
precon_thblock integer in ブロックヤコビ前処理,加法シュワルツ前処理のブロック数
independ_nvec integer in 細分化前処理で独立に計算する対角ブロックのサイズ
nBlock integer in サイクリックにスレッドに割り当てる行数
iret integer out エラーコード(0:正常終了)


6.8 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, precon_thblock, independ_nvec, nBlock, iret )


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

引数名(次元) 入力/出力 説明
icomm integer in MPIコミュニケータ
vecx(n) double precision in/out in:反復の初期値,out:解ベクトル
vecb(n) double precision in 右辺ベクトル
n integer*4 / integer*8 in 各プロセスが担当するベクトルの大きさ
gn integer*4 / integer*8 in 全体のベクトルの大きさ
istart integer*4 / integer*8 in 各プロセスが担当する行列の開始行
diaA(n*nnd) double precision in DIA形式に格納された行列の非零要素
offset(nnd) integer*4 / integer*8 in DIA形式に格納された行列の対角からのオフセット
nnd integer*4 / integer*8 in DIA形式に格納された行列の対角要素の本数
ipreflag integer in 前処理指定フラグ(0:無し,1:ポイントヤコビ,2:ブロックヤコビ,3:加法シュワルツ)
ilu_method integer in 不完全LU分解指定フラグ(0:ILU(0),1:D-ILU(対角成分が一致),2:D-ILU(行内の要素和が一致),加法シュワルツは1,2のみ利用可能)
addL integer in 加法シュワルツの重なり幅
iflagAS integer in 加法シュワルツ指定フラグ(1:BASIC,2:RESTRICT,3:INTERPOLATE,4:NONE)
itrmax integer in/out in:最大反復回数,out:反復回数
rtolmax double precision in 収束判定条件(相対残差ノルム)
reshistory(itrmax) double precision out 残差履歴
iovlflag integer in 通信隠蔽フラグ(0:無し,1:全プロセス,2:各プロセス)
precon_thblock integer in ブロックヤコビ前処理,加法シュワルツ前処理のブロック数
independ_nvec integer in 細分化前処理で独立に計算する対角ブロックのサイズ
nBlock integer in サイクリックにスレッドに割り当てる行数
iret integer out エラーコード(0:正常終了)


6.9 parcel_dgmres_dia


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


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

引数名(次元) 入力/出力 説明
icomm integer in MPIコミュニケータ
vecx(n) double precision in/out in:反復の初期値,out:解ベクトル
vecb(n) double precision in 右辺ベクトル
n integer*4 / integer*8 in 各プロセスが担当するベクトルの大きさ
gn integer*4 / integer*8 in 全体のベクトルの大きさ
istart integer*4 / integer*8 in 各プロセスが担当する行列の開始行
diaA(n*nnd) double precision in DIA形式に格納された行列の非零要素
offset(nnd) integer*4 / integer*8 in DIA形式に格納された行列の対角からのオフセット
nnd integer*4 / integer*8 in DIA形式に格納された行列の対角要素の本数
ipreflag integer in 前処理指定フラグ(0:無し,1:ポイントヤコビ,2:ブロックヤコビ,3:加法シュワルツ)
ilu_method integer in 不完全LU分解指定フラグ(0:ILU(0),1:D-ILU(対角成分が一致),2:D-ILU(行内の要素和が一致),加法シュワルツは1,2のみ利用可能)
addL integer in 加法シュワルツの重なり幅
iflagAS integer in 加法シュワルツ指定フラグ(1:BASIC,2:RESTRICT,3:INTERPOLATE,4:NONE)
itrmax integer in/out in:最大反復回数,out:反復回数
rtolmax double precision in 収束判定条件(相対残差ノルム)
reshistory(itrmax) double precision out 残差履歴
iovlflag integer in 通信隠蔽フラグ(0:無し,1:全プロセス,2:各プロセス)
precon_thblock integer in ブロックヤコビ前処理,加法シュワルツ前処理のブロック数
independ_nvec integer in 細分化前処理で独立に計算する対角ブロックのサイズ
nBlock integer in サイクリックにスレッドに割り当てる行数
gmres_m integer in リスタートするまでの反復回数
gmres_GSflag integer in 直交化フラグ(1:MGS,2:CGS)
iret integer out エラーコード(0:正常終了)


6.10 parcel_dcagmres_dia


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


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

引数名(次元) 入力/出力 説明
icomm integer in MPIコミュニケータ
vecx(n) double precision in/out in:反復の初期値,out:解ベクトル
vecb(n) double precision in 右辺ベクトル
n integer*4 / integer*8 in 各プロセスが担当するベクトルの大きさ
gn integer*4 / integer*8 in 全体のベクトルの大きさ
istart integer*4 / integer*8 in 各プロセスが担当する行列の開始行
diaA(n*nnd) double precision in DIA形式に格納された行列の非零要素
offset(nnd) integer*4 / integer*8 in DIA形式に格納された行列の対角からのオフセット
nnd integer*4 / integer*8 in DIA形式に格納された行列の対角要素の本数
ipreflag integer in 前処理指定フラグ(0:無し,1:ポイントヤコビ,2:ブロックヤコビ,3:加法シュワルツ)
ilu_method integer in 不完全LU分解指定フラグ(0:ILU(0),1:D-ILU(対角成分が一致),2:D-ILU(行内の要素和が一致),加法シュワルツは1,2のみ利用可能)
addL integer in 加法シュワルツの重なり幅
iflagAS integer in 加法シュワルツ指定フラグ(1:BASIC,2:RESTRICT,3:INTERPOLATE,4:NONE)
itrmax integer in/out in:最大反復回数,out:反復回数
rtolmax double precision in 収束判定条件(相対残差ノルム)
reshistory(itrmax) double precision out 残差履歴
iovlflag integer in 通信隠蔽フラグ(0:無し,1:全プロセス,2:各プロセス)
precon_thblock integer in ブロックヤコビ前処理,加法シュワルツ前処理のブロック数
independ_nvec integer in 細分化前処理で独立に計算する対角ブロックのサイズ
nBlock integer in サイクリックにスレッドに割り当てる行数
cagmres_sstep integer in CA-GMRES法の省通信ステップ
cagmres_tstep integer in CA-GMRES法のリスタートステップ
cagmres_basis integer in CA-GMRES法の基底生成フラグ(0:単基底,1:ニュートン基底)
cagmres_QRflag integer in CA-GMRES法のQR分解フラグ(1:MGS,2:CGS,3:TSQR,4:CholeskyQR,5:CholeskyQR2)
iret integer out エラーコード(0:正常終了)


6.11 parcel_dcbcg_dia


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


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

引数名(次元) 入力/出力 説明
icomm integer in MPIコミュニケータ
vecx(n) double precision in/out in:反復の初期値,out:解ベクトル
vecb(n) double precision in 右辺ベクトル
n integer*4 / integer*8 in 各プロセスが担当するベクトルの大きさ
gn integer*4 / integer*8 in 全体のベクトルの大きさ
istart integer*4 / integer*8 in 各プロセスが担当する行列の開始行
diaA(n*nnd) double precision in DIA形式に格納された行列の非零要素
offset(nnd) integer*4 / integer*8 in DIA形式に格納された行列の対角からのオフセット
nnd integer*4 / integer*8 in DIA形式に格納された行列の対角要素の本数
ipreflag integer in 前処理指定フラグ(0:無し,1:ポイントヤコビ,2:ブロックヤコビ,3:加法シュワルツ)
ilu_method integer in 不完全LU分解指定フラグ(0:ILU(0),1:D-ILU(対角成分が一致),2:D-ILU(行内の要素和が一致),加法シュワルツは1,2のみ利用可能)
addL integer in 加法シュワルツの重なり幅
iflagAS integer in 加法シュワルツ指定フラグ(1:BASIC,2:RESTRICT,3:INTERPOLATE,4:NONE)
itrmax integer in/out in:最大反復回数,out:反復回数
rtolmax double precision in 収束判定条件(相対残差ノルム)
reshistory(itrmax) double precision out 残差履歴
iovlflag integer in 通信隠蔽フラグ(0:無し,1:全プロセス,2:各プロセス)
precon_thblock integer in ブロックヤコビ前処理,加法シュワルツ前処理のブロック数
independ_nvec integer in 細分化前処理で独立に計算する対角ブロックのサイズ
nBlock integer in サイクリックにスレッドに割り当てる行数
cbcg_kstep integer in CBCG法の省通信ステップ
cbcg_Eigenflag integer in 固有値計算フラグ(1:べき乗法,2:CA-ARNOLDI)
power_method_itrmax integer in べき乗法の最大反復回数
caarnoldi_sstep integer in CA-ARNOLDI法の省通信ステップ
caarnoldi_tstep integer in CA-ARNOLDI法のリスタートステップ
caarnoldi_basis integer in CA-ARNOLDI法の基底生成フラグ(0:単基底,1:ニュートン基底)
caarnoldi_QRflag integer in CA-ARNOLDI法のQR分解フラグ(1:MGS,2:CGS,3:TSQR,4:CholeskyQR,5:CholeskyQR2)
iret integer out エラーコード(0:正常終了)


6.12 parcel_dcaarnoldi_dia


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形式で格納.

引数名(次元) 入力/出力 説明
icomm integer in MPIコミュニケータ
vecx(n) double precision in/out in:反復の初期値,out:解ベクトル
vecb(n) double precision in 右辺ベクトル
n integer*4 / integer*8 in 各プロセスが担当するベクトルの大きさ
gn integer*4 / integer*8 in 全体のベクトルの大きさ
istart integer*4 / integer*8 in 各プロセスが担当する行列の開始行
diaA(n*nnd) double precision in DIA形式に格納された行列の非零要素
offset(nnd) integer*4 / integer*8 in DIA形式に格納された行列の対角からのオフセット
nnd integer*4 / integer*8 in DIA形式に格納された行列の対角要素の本数
ipreflag integer in 前処理指定フラグ(0:無し,1:ポイントヤコビ,2:ブロックヤコビ,3:加法シュワルツ)
ilu_method integer in 不完全LU分解指定フラグ(0:ILU(0),1:D-ILU(対角成分が一致),2:D-ILU(行内の要素和が一致),加法シュワルツは1,2のみ利用可能)
addL integer in 加法シュワルツの重なり幅
iflagAS integer in 加法シュワルツ指定フラグ(1:BASIC,2:RESTRICT,3:INTERPOLATE,4:NONE)
itrmax integer in/out in:最大反復回数,out:反復回数
iovlflag integer in 通信隠蔽フラグ(0:無し,1:全プロセス,2:各プロセス)
precon_thblock integer in ブロックヤコビ前処理,加法シュワルツ前処理のブロック数
independ_nvec integer in 細分化前処理で独立に計算する対角ブロックのサイズ
nBlock integer in サイクリックにスレッドに割り当てる行数
caarnoldi_sstep integer in CA-ARNOLDI法の省通信ステップ数
caarnoldi_tstep integer in CA-ARNOLDI法のリスタートステップ数
caarnoldi_basis integer in CA-ARNOLDI法の基底生成フラグ(0:単基底,1:ニュートン基底)
caarnoldi_QRflag integer in CA-ARNOLDI法のQR分解フラグ(1:MGS,2:CGS,3:TSQR,4:CholeskyQR,5:CholeskyQR2)
Evalr(caarnoldi_sstep*caarnoldi_tstep) double precision out 固有値の実数部
Evali(caarnoldi_sstep*caarnoldi_tstep) double precision out 固有値の虚数部
Evec(n*caarnoldi_sstep*caarnoldi_tstep) double precision out 固有ベクトル
Eerr(caarnoldi_sstep*caarnoldi_tstep) double precision out 固有値及び固有ベクトルのエラーノルム \({\left\| Ax - \lambda x \right\|}/{\left\| \text{λx} \right\|}\)
EmaxID integer out 最大固有値のID番号
EminID integer out 最小固有値のID番号
iret integer out エラーコード(0:正常終了)


6.13 parcel_dcg_ddm


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, independ_nvec, nBlock, iret )


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

引数名(次元) 入力/出力 説明
icomm integer in MPIコミュニケータ
x(m) double precision in/out in:反復の初期値,out:解ベクトル
b(m) double precision in 右辺ベクトル
m integer in 各プロセスが担当するベクトルの大きさ
nnd integer in プロセスが担当する係数行列要素の非零成分をもつ対角の本数
ioff_dia(nnd) integer in DDM形式に格納された行列の対角からのオフセット
val_dia(m*nnd) double precision in DDM形式に格納された行列の非ゼロ要素
num_neighbor_ranks integer in 隣接するプロセス数
neighbor_ranks(num_neighbor_ranks) integer in 隣接するプロセス番号のリスト
margin integer in ioff_gridの絶対値の最大値
ncomp_grids integer in 構造格子の次元数
ndiv_grids(ncomp_grids) integer in 構造格子の各次元の分割数
num_grids(ncomp_grids) integer in プロセスが担当する格子サイズ
ioff_grids(nnd,ncomp_grids) integer in ステンシルの各次元方向の成分
div_direc_th integer in スレッド並列における領域分割の方向(1:x方向,2:y方向,3:z方向)
ipreflag integer in 前処理指定フラグ(0:無し,1:ポイントヤコビ,2:ブロックヤコビ,3:加法シュワルツ
addL integer in 加法シュワルツの重なり幅
iflagAS integer in 加法シュワルツ指定フラグ(1:BASIC,2:RESTRICT,3:INTERPOLATE,4:NONE)
itrmax integer in/out in:最大反復回数,out:収束回数
rtolmax double precision in 収束判定条件(相対残差ノルム)
reshistory(itrmax) double precision out 残差履歴
iovlflag integer in 通信隠蔽フラグ(0:無し,1:全プロセス,2:各プロセス)
precon_thblock integer in ブロックヤコビ前処理,加法シュワルツ前処理のブロック数
independ_nvec integer in 細分化前処理で独立に計算する対角ブロックのサイズ
nBlock integer in サイクリックにスレッドに割り当てる行数
iret integer out エラーコード(0:正常終了)


6.14 parcel_dbicgstab_ddm


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, independ_nvec, nBlock, iret )


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

引数名(次元) 入力/出力 説明
icomm integer in MPIコミュニケータ
x(m) double precision in/out in:反復の初期値,out:解ベクトル
b(m) double precision in 右辺ベクトル
m integer in 各プロセスが担当するベクトルの大きさ
nnd integer in プロセスが担当する係数行列要素の非零成分をもつ対角の本数
ioff_dia(nnd) integer in DDM形式に格納された行列の対角からのオフセット
val_dia(m*nnd) double precision in DDM形式に格納された行列の非ゼロ要素
num_neighbor_ranks integer in 隣接するプロセス数
neighbor_ranks(num_neighbor_ranks) integer in 隣接するプロセス番号のリスト
margin integer in ioff_gridの絶対値の最大値
ncomp_grids integer in 構造格子の次元数
ndiv_grids(ncomp_grids) integer in 構造格子の各次元の分割数
num_grids(ncomp_grids) integer in プロセスが担当する格子サイズ
ioff_grids(nnd,ncomp_grids) integer in ステンシルの各次元方向の成分
div_direc_th integer in スレッド並列における領域分割の方向(1:x方向,2:y方向,3:z方向)
ipreflag integer in 前処理指定フラグ(0:無し,1:ポイントヤコビ,2:ブロックヤコビ,3:加法シュワルツ
addL integer in 加法シュワルツの重なり幅
iflagAS integer in 加法シュワルツ指定フラグ(1:BASIC,2:RESTRICT,3:INTERPOLATE,4:NONE)
itrmax integer in/out in:最大反復回数,out:収束回数
rtolmax double precision in 収束判定条件(相対残差ノルム)
reshistory(itrmax) double precision out 残差履歴
iovlflag integer in 通信隠蔽フラグ(0:無し,1:全プロセス,2:各プロセス)
precon_thblock integer in ブロックヤコビ前処理,加法シュワルツ前処理のブロック数
independ_nvec integer in 細分化前処理で独立に計算する対角ブロックのサイズ
nBlock integer in サイクリックにスレッドに割り当てる行数
iret integer out エラーコード(0:正常終了)


6.15 parcel_dgmres_ddm


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, independ_nvec, nBlock, gmres_m, gmres_GSflag, iret )


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

引数名(次元) 入力/出力 説明
icomm integer in MPIコミュニケータ
x(m) double precision in/out in:反復の初期値,out:解ベクトル
b(m) double precision in 右辺ベクトル
m integer in 各プロセスが担当するベクトルの大きさ
nnd integer in プロセスが担当する係数行列要素の非零成分をもつ対角の本数
ioff_dia(nnd) integer in DDM形式に格納された行列の対角からのオフセット
val_dia(m*nnd) double precision in DDM形式に格納された行列の非ゼロ要素
num_neighbor_ranks integer in 隣接するプロセス数
neighbor_ranks(num_neighbor_ranks) integer in 隣接するプロセス番号のリスト
margin integer in ioff_gridの絶対値の最大値
ncomp_grids integer in 構造格子の次元数
ndiv_grids(ncomp_grids) integer in 構造格子の各次元の分割数
num_grids(ncomp_grids) integer in プロセスが担当する格子サイズ
ioff_grids(nnd,ncomp_grids) integer in ステンシルの各次元方向の成分
div_direc_th integer in スレッド並列における領域分割の方向(1:x方向,2:y方向,3:z方向)
ipreflag integer in 前処理指定フラグ(0:無し,1:ポイントヤコビ,2:ブロックヤコビ,3:加法シュワルツ
addL integer in 加法シュワルツの重なり幅
iflagAS integer in 加法シュワルツ指定フラグ(1:BASIC,2:RESTRICT,3:INTERPOLATE,4:NONE)
itrmax integer in/out in:最大反復回数,out:収束回数
rtolmax double precision in 収束判定条件(相対残差ノルム)
reshistory(itrmax) double precision out 残差履歴
iovlflag integer in 通信隠蔽フラグ(0:無し,1:全プロセス,2:各プロセス)
precon_thblock integer in ブロックヤコビ前処理,加法シュワルツ前処理のブロック数
independ_nvec integer in 細分化前処理で独立に計算する対角ブロックのサイズ
nBlock integer in サイクリックにスレッドに割り当てる行数
gmres_m integer in リスタートするまでの反復回数
gmres_GSflag integer in 直交化フラグ(1:MGS,2:CGS)
iret integer out エラーコード(0:正常終了)


6.16 parcel_dcagmres_ddm


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, independ_nvec, nBlock, cagmres_sstep, cagmres_tstep, cagmres_basis, cagmres_QRflag, iret )


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

引数名(次元) 入力/出力 説明
icomm integer in MPIコミュニケータ
x(m) double precision in/out in:反復の初期値,out:解ベクトル
b(m) double precision in 右辺ベクトル
m integer in 各プロセスが担当するベクトルの大きさ
nnd integer in プロセスが担当する係数行列要素の非零成分をもつ対角の本数
ioff_dia(nnd) integer in DDM形式に格納された行列の対角からのオフセット
val_dia(m*nnd) double precision in DDM形式に格納された行列の非ゼロ要素
num_neighbor_ranks integer in 隣接するプロセス数
neighbor_ranks(num_neighbor_ranks) integer in 隣接するプロセス番号のリスト
margin integer in ioff_gridの絶対値の最大値
ncomp_grids integer in 構造格子の次元数
ndiv_grids(ncomp_grids) integer in 構造格子の各次元の分割数
num_grids(ncomp_grids) integer in プロセスが担当する格子サイズ
ioff_grids(nnd,ncomp_grids) integer in ステンシルの各次元方向の成分
div_direc_th integer in スレッド並列における領域分割の方向(1:x方向,2:y方向,3:z方向)
ipreflag integer in 前処理指定フラグ(0:無し,1:ポイントヤコビ,2:ブロックヤコビ,3:加法シュワルツ
addL integer in 加法シュワルツの重なり幅
iflagAS integer in 加法シュワルツ指定フラグ(1:BASIC,2:RESTRICT,3:INTERPOLATE,4:NONE)
itrmax integer in/out in:最大反復回数,out:収束回数
rtolmax double precision in 収束判定条件(相対残差ノルム)
reshistory(itrmax) double precision out 残差履歴
iovlflag integer in 通信隠蔽フラグ(0:無し,1:全プロセス,2:各プロセス)
precon_thblock integer in ブロックヤコビ前処理,加法シュワルツ前処理のブロック数
independ_nvec integer in 細分化前処理で独立に計算する対角ブロックのサイズ
nBlock integer in サイクリックにスレッドに割り当てる行数
cagmres_sstep integer in CA-GMRES法の省通信ステップ数
cagmres_tstep integer in CA-GMRES法のリスタートステップ数
cagmres_basis integer in CA-GMRES法の基底生成フラグ(0:単基底,1:ニュートン基底)
cagmres_QRflag integer in CA-GMRES法のQR分解フラグ(1:MGS,2:CGS,3:TSQR,4:CholeskyQR,5:CholeskyQR2)
iret integer out エラーコード(0:正常終了)


6.17 parcel_dcbcg_ddm


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, independ_nvec, nBlock, cbcg_kstep, cbcg_Eigenflag, power_method_itrmax, caarnoldi_sstep, caarnoldi_tstep, caarnoldi_basis, caarnoldi_QRflag, iret )


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

引数名(次元) 入力/出力 説明
icomm integer in MPIコミュニケータ
x(m) double precision in/out in:反復の初期値,out:解ベクトル
b(m) double precision in 右辺ベクトル
m integer in 各プロセスが担当するベクトルの大きさ
nnd integer in プロセスが担当する係数行列要素の非零成分をもつ対角の本数
ioff_dia(nnd) integer in DDM形式に格納された行列の対角からのオフセット
val_dia(m*nnd) double precision in DDM形式に格納された行列の非ゼロ要素
num_neighbor_ranks integer in 隣接するプロセス数
neighbor_ranks(num_neighbor_ranks) integer in 隣接するプロセス番号のリスト
margin integer in ioff_gridの絶対値の最大値
ncomp_grids integer in 構造格子の次元数
ndiv_grids(ncomp_grids) integer in 構造格子の各次元の分割数
num_grids(ncomp_grids) integer in プロセスが担当する格子サイズ
ioff_grids(nnd,ncomp_grids) integer in ステンシルの各次元方向の成分
div_direc_th integer in スレッド並列における領域分割の方向(1:x方向,2:y方向,3:z方向)
ipreflag integer in 前処理指定フラグ(0:無し,1:ポイントヤコビ,2:ブロックヤコビ,3:加法シュワルツ
addL integer in 加法シュワルツの重なり幅
iflagAS integer in 加法シュワルツ指定フラグ(1:BASIC,2:RESTRICT,3:INTERPOLATE,4:NONE)
itrmax integer in/out in:最大反復回数,out:収束回数
rtolmax double precision in 収束判定条件(相対残差ノルム)
reshistory(itrmax) double precision out 残差履歴
iovlflag integer in 通信隠蔽フラグ(0:無し,1:全プロセス,2:各プロセス)
precon_thblock integer in ブロックヤコビ前処理,加法シュワルツ前処理のブロック数
independ_nvec integer in 細分化前処理で独立に計算する対角ブロックのサイズ
nBlock integer in サイクリックにスレッドに割り当てる行数
cbcg_kstep integer in CBCG法の省通信ステップ
cbcg_Eigenflag integer in 固有値計算フラグ(1:べき乗法,2:CA-ARNOLDI)
power_method_itrmax integer in べき乗法の反復回数
caarnoldi_sstep integer in CA-ARNOLDI法の省通信ステップ数
caarnoldi_tstep integer in CA-ARNOLDI法のリスタートステップ数
caarnoldi_basis integer in CA-ARNOLDI法の基底生成フラグ(0:単基底,1:ニュートン基底)
caarnoldi_QRflag integer in CA-ARNOLDI法のQR分解フラグ(1:MGS,2:CGS,3:TSQR,4:CholeskyQR,5:CholeskyQR2)
iret integer out エラーコード(0:正常終了)


6.18 parcel_dcaarnoldi_ddm


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形式で格納.

引数名(次元) 入力/出力 説明
icomm integer in MPIコミュニケータ
x(m) double precision in/out in:反復の初期値,out:解ベクトル
b(m) double precision in 右辺ベクトル
m integer in 各プロセスが担当するベクトルの大きさ
nnd integer in プロセスが担当する係数行列要素の非零成分をもつ対角の本数
ioff_dia(nnd) integer in DDM形式に格納された行列の対角からのオフセット
val_dia(m*nnd) double precision in DDM形式に格納された行列の非ゼロ要素
num_neighbor_ranks integer in 隣接するプロセス数
neighbor_ranks(num_neighbor_ranks) integer in 隣接するプロセス番号のリスト
margin integer in ioff_gridの絶対値の最大値
ncomp_grids integer in 構造格子の次元数
ndiv_grids(ncomp_grids) integer in 構造格子の各次元の分割数
num_grids(ncomp_grids) integer in プロセスが担当する格子サイズ
ioff_grids(nnd,ncomp_grids) integer in ステンシルの各次元方向の成分
div_direc_th integer in スレッド並列における領域分割の方向(1:x方向,2:y方向,3:z方向)
ipreflag integer in 前処理指定フラグ(0:無し,1:ポイントヤコビ,2:ブロックヤコビ,3:加法シュワルツ
addL integer in 加法シュワルツの重なり幅
iflagAS integer in 加法シュワルツ指定フラグ(1:BASIC,2:RESTRICT,3:INTERPOLATE,4:NONE)
iovlflag integer in 通信隠蔽フラグ(0:無し,1:全プロセス,2:各プロセス)
precon_thblock integer in ブロックヤコビ前処理,加法シュワルツ前処理のブロック数
independ_nvec integer in 細分化前処理で独立に計算する対角ブロックのサイズ
nBlock integer in サイクリックにスレッドに割り当てる行数
caarnoldi_sstep integer in CA-ARNOLDI法の省通信ステップ数
caarnoldi_tstep integer in CA-ARNOLDI法のリスタートステップ数
caarnoldi_basis integer in CA-ARNOLDI法の基底生成フラグ(0:単基底,1:ニュートン基底)
caarnoldi_QRflag integer in CA-ARNOLDI法のQR分解フラグ(1:MGS,2:CGS,3:TSQR,4:CholeskyQR,5:CholeskyQR2)
Evalr(caarnoldi_sstep*caarnoldi_tstep) double precision out 固有値の実数部
Evali(caarnoldi_sstep*caarnoldi_tstep) double precision out 固有値の虚数部
Evec(m*caarnoldi_sstep*caarnoldi_tstep) double precision out 固有ベクトル
Eerr(caarnoldi_sstep*caarnoldi_tstep) double precision out 固有値及び固有ベクトルのエラーノルム \({\left\| Ax - \lambda x \right\|}/{\left\| \text{λx} \right\|}\)
EmaxID integer out 最大固有値のID番号
EminID integer out 最小固有値のID番号
iret integer out エラーコード(0:正常終了)


6.19 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, precon_thblock, independ_nvec, nBlock, precision_A, precision_b, precision_x, precision_precon, iret )


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

引数名(次元) 入力/出力 説明
icomm integer in MPIコミュニケータ
vecx_hi(n) double precision in/out in:反復の初期値(上位ビット),out:解ベクトル(上位ビット)
vecx_lo(n) double precision in/out in:反復の初期値(下位ビット),out:解ベクトル(下位ビット)
vecb_hi(n) double precision in 右辺ベクトル(上位ビット)
vecb_lo(n) double precision in 右辺ベクトル(下位ビット)
n integer*4 / integer*8 in 各プロセスが担当するベクトルの大きさ
gn integer*4 / integer*8 in 全体のベクトルの大きさ
nnz integer*4 / integer*8 in 各プロセスが担当する係数行列の非零要素数
istart integer*4 / integer*8 in 各プロセスが担当する行列の開始行
crsA_hi(nnz) double precision in CRS 形式に格納された行列の非零要素(上位ビット)
crsA_lo(nnz) double precision in CRS 形式に格納された行列の非零要素(下位ビット)
crsRow_ptr(n+1) integer*4 / integer*8 in CRS 形式のポインタテーブル
crsCol(nnz) integer*4 / integer*8 in CRS 形式に格納された分割行列の非零要素の列番号
ipreflag integer in 前処理指定フラグ(0:無し,1:ポイントヤコビ,2:ブロックヤコビ,3:加法シュワルツ)
ilu_method integer in 不完全LU分解指定フラグ(0:ILU(0),1:D-ILU(対角成分が一致),2:D-ILU(行内の要素和が一致),加法シュワルツは1,2のみ利用可能)
addL integer in 加法シュワルツの重なり幅
iflagAS integer in 加法シュワルツ指定フラグ(1:BASIC,2:RESTRICT,3:INTERPOLATE,4:NONE)
itrmax integer in/out in:最大反復回数,out:反復回数
rtolmax double precision in 収束判定条件(相対残差ノルム)
reshistory(itrmax) double precision out 残差履歴
iovlflag integer in 通信隠蔽フラグ(0:無し,1:全プロセス,2:各プロセス)
precon_thblock integer in ブロックヤコビ前処理,加法シュワルツ前処理のブロック数
independ_nvec integer in 細分化前処理で独立に計算する対角ブロックのサイズ
nBlock integer in サイクリックにスレッドに割り当てる行数
precision_A integer in 行列の非零要素の精度(1:倍精度,2:4倍精度)
precision_b integer in 右辺ベクトルの精度(1:倍精度,2:4倍精度)
precision_x integer in 反復の初期値,解ベクトルの精度(1:倍精度,2:4倍精度)
precision_precon integer in 前処理行列の精度(1:倍精度,2:4倍精度)
iret integer out エラーコード(0:正常終了)


6.20 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, precon_thblock, independ_nvec, nBlock, precision_A, precision_b, precision_x, precision_precon, iret )


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

引数名(次元) 入力/出力 説明
icomm integer in MPIコミュニケータ
vecx_hi(n) double precision in/out in:反復の初期値(上位ビット),out:解ベクトル(上位ビット)
vecx_lo(n) double precision in/out in:反復の初期値(下位ビット),out:解ベクトル(下位ビット)
vecb_hi(n) double precision in 右辺ベクトル(上位ビット)
vecb_lo(n) double precision in 右辺ベクトル(下位ビット)
n integer*4 / integer*8 in 各プロセスが担当するベクトルの大きさ
gn integer*4 / integer*8 in 全体のベクトルの大きさ
nnz integer*4 / integer*8 in 各プロセスが担当する係数行列の非零要素数
istart integer*4 / integer*8 in 各プロセスが担当する行列の開始行
crsA_hi(nnz) double precision in CRS 形式に格納された行列の非零要素(上位ビット)
crsA_lo(nnz) double precision in CRS 形式に格納された行列の非零要素(下位ビット)
crsRow_ptr(n+1) integer*4 / integer*8 in CRS 形式のポインタテーブル
crsCol(nnz) integer*4 / integer*8 in CRS 形式に格納された分割行列の非零要素の列番号
ipreflag integer in 前処理指定フラグ(0:無し,1:ポイントヤコビ,2:ブロックヤコビ,3:加法シュワルツ)
ilu_method integer in 不完全LU分解指定フラグ(0:ILU(0),1:D-ILU(対角成分が一致),2:D-ILU(行内の要素和が一致),加法シュワルツは1,2のみ利用可能)
addL integer in 加法シュワルツの重なり幅
iflagAS integer in 加法シュワルツ指定フラグ(1:BASIC,2:RESTRICT,3:INTERPOLATE,4:NONE)
itrmax integer in/out in:最大反復回数,out:反復回数
rtolmax double precision in 収束判定条件(相対残差ノルム)
reshistory(itrmax) double precision out 残差履歴
iovlflag integer in 通信隠蔽フラグ(0:無し,1:全プロセス,2:各プロセス)
precon_thblock integer in ブロックヤコビ前処理,加法シュワルツ前処理のブロック数
independ_nvec integer in 細分化前処理で独立に計算する対角ブロックのサイズ
nBlock integer in サイクリックにスレッドに割り当てる行数
precision_A integer in 行列の非零要素の精度(1:倍精度,2:4倍精度)
precision_b integer in 右辺ベクトルの精度(1:倍精度,2:4倍精度)
precision_x integer in 反復の初期値,解ベクトルの精度(1:倍精度,2:4倍精度)
precision_precon integer in 前処理行列の精度(1:倍精度,2:4倍精度)
iret integer out エラーコード(0:正常終了)


6.21 parcel_ddgmres


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形式で格納.

引数名(次元) 入力/出力 説明
icomm integer in MPIコミュニケータ
vecx_hi(n) double precision in/out in:反復の初期値(上位ビット),out:解ベクトル(上位ビット)
vecx_lo(n) double precision in/out in:反復の初期値(下位ビット),out:解ベクトル(下位ビット)
vecb_hi(n) double precision in 右辺ベクトル(上位ビット)
vecb_lo(n) double precision in 右辺ベクトル(下位ビット)
n integer*4 / integer*8 in 各プロセスが担当するベクトルの大きさ
gn integer*4 / integer*8 in 全体のベクトルの大きさ
nnz integer*4 / integer*8 in 各プロセスが担当する係数行列の非零要素数
istart integer*4 / integer*8 in 各プロセスが担当する行列の開始行
crsA_hi(nnz) double precision in CRS 形式に格納された行列の非零要素(上位ビット)
crsA_lo(nnz) double precision in CRS 形式に格納された行列の非零要素(下位ビット)
crsRow_ptr(n+1) integer*4 / integer*8 in CRS 形式のポインタテーブル
crsCol(nnz) integer*4 / integer*8 in CRS 形式に格納された分割行列の非零要素の列番号
ipreflag integer in 前処理指定フラグ(0:無し,1:ポイントヤコビ,2:ブロックヤコビ,3:加法シュワルツ)
ilu_method integer in 不完全LU分解指定フラグ(0:ILU(0),1:D-ILU(対角成分が一致),2:D-ILU(行内の要素和が一致),加法シュワルツは1,2のみ利用可能)
addL integer in 加法シュワルツの重なり幅
iflagAS integer in 加法シュワルツ指定フラグ(1:BASIC,2:RESTRICT,3:INTERPOLATE,4:NONE)
itrmax integer in/out in:最大反復回数,out:反復回数
rtolmax double precision in 収束判定条件(相対残差ノルム)
reshistory(itrmax) double precision out 残差履歴
iovlflag integer in 通信隠蔽フラグ(0:無し,1:全プロセス,2:各プロセス)
precon_thblock integer in ブロックヤコビ前処理,加法シュワルツ前処理のブロック数
independ_nvec integer in 細分化前処理で独立に計算する対角ブロックのサイズ
nBlock integer in サイクリックにスレッドに割り当てる行数
gmres_m integer in リスタートするまでの反復回数
gmres_GSflag integer in 直交化フラグ(1:MGS,2:CGS)
precision_A integer in 行列の非零要素の精度(1:倍精度,2:4倍精度)
precision_b integer in 右辺ベクトルの精度(1:倍精度,2:4倍精度)
precision_x integer in 反復の初期値,解ベクトルの精度(1:倍精度,2:4倍精度)
precision_precon integer in 前処理行列の精度(1:倍精度,2:4倍精度)
iret integer out エラーコード(0:正常終了)


6.22 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, precon_thblock, independ_nvec, nBlock, precision_A, precision_b, precision_x, precision_precon, iret )


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

引数名(次元) 入力/出力 説明
icomm integer in MPIコミュニケータ
vecx_hi(n) double precision in/out in:反復の初期値(上位ビット),out:解ベクトル(上位ビット)
vecx_lo(n) double precision in/out in:反復の初期値(下位ビット),out:解ベクトル(下位ビット)
vecb_hi(n) double precision in 右辺ベクトル(上位ビット)
vecb_lo(n) double precision in 右辺ベクトル(下位ビット)
n integer*4 / integer*8 in 各プロセスが担当するベクトルの大きさ
gn integer*4 / integer*8 in 全体のベクトルの大きさ
istart integer*4 / integer*8 in 各プロセスが担当する行列の開始行
diaA_hi(n*nnd) double precision in DIA形式に格納された行列の非零要素(上位ビット)
diaA_lo(n*nnd) double precision in DIA形式に格納された行列の非零要素(下位ビット)
offset(nnd) integer*4 / integer*8 in DIA形式に格納された行列の対角からのオフセット
nnd integer*4 / integer*8 in DIA形式に格納された行列の対角要素の本数
ipreflag integer in 前処理指定フラグ(0:無し,1:ポイントヤコビ,2:ブロックヤコビ,3:加法シュワルツ)
ilu_method integer in 不完全LU分解指定フラグ(0:ILU(0),1:D-ILU(対角成分が一致),2:D-ILU(行内の要素和が一致),加法シュワルツは1,2のみ利用可能)
addL integer in 加法シュワルツの重なり幅
iflagAS integer in 加法シュワルツ指定フラグ(1:BASIC,2:RESTRICT,3:INTERPOLATE,4:NONE)
itrmax integer in/out in:最大反復回数,out:反復回数
rtolmax double precision in 収束判定条件(相対残差ノルム)
reshistory(itrmax) double precision out 残差履歴
iovlflag integer in 通信隠蔽フラグ(0:無し,1:全プロセス,2:各プロセス)
precon_thblock integer in ブロックヤコビ前処理,加法シュワルツ前処理のブロック数
independ_nvec integer in 細分化前処理で独立に計算する対角ブロックのサイズ
nBlock integer in サイクリックにスレッドに割り当てる行数
precision_A integer in 行列の非零要素の精度(1:倍精度,2:4倍精度)
precision_b integer in 右辺ベクトルの精度(1:倍精度,2:4倍精度)
precision_x integer in 反復の初期値,解ベクトルの精度(1:倍精度,2:4倍精度)
precision_precon integer in 前処理行列の精度(1:倍精度,2:4倍精度)
iret integer out エラーコード(0:正常終了)


6.23 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, precon_thblock, independ_nvec, nBlock, precision_A, precision_b, precision_x, precision_precon, iret )


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

引数名(次元) 入力/出力 説明
icomm integer in MPIコミュニケータ
vecx_hi(n) double precision in/out in:反復の初期値(上位ビット),out:解ベクトル(上位ビット)
vecx_lo(n) double precision in/out in:反復の初期値(下位ビット),out:解ベクトル(下位ビット)
vecb_hi(n) double precision in 右辺ベクトル(上位ビット)
vecb_lo(n) double precision in 右辺ベクトル(下位ビット)
n integer*4 / integer*8 in 各プロセスが担当するベクトルの大きさ
gn integer*4 / integer*8 in 全体のベクトルの大きさ
istart integer*4 / integer*8 in 各プロセスが担当する行列の開始行
diaA_hi(n*nnd) double precision in DIA形式に格納された行列の非零要素(上位ビット)
diaA_lo(n*nnd) double precision in DIA形式に格納された行列の非零要素(下位ビット)
offset(nnd) integer*4 / integer*8 in DIA形式に格納された行列の対角からのオフセット
nnd integer*4 / integer*8 in DIA形式に格納された行列の対角要素の本数
ipreflag integer in 前処理指定フラグ(0:無し,1:ポイントヤコビ,2:ブロックヤコビ,3:加法シュワルツ)
ilu_method integer in 不完全LU分解指定フラグ(0:ILU(0),1:D-ILU(対角成分が一致),2:D-ILU(行内の要素和が一致),加法シュワルツは1,2のみ利用可能)
addL integer in 加法シュワルツの重なり幅
iflagAS integer in 加法シュワルツ指定フラグ(1:BASIC,2:RESTRICT,3:INTERPOLATE,4:NONE)
itrmax integer in/out in:最大反復回数,out:反復回数
rtolmax double precision in 収束判定条件(相対残差ノルム)
reshistory(itrmax) double precision out 残差履歴
iovlflag integer in 通信隠蔽フラグ(0:無し,1:全プロセス,2:各プロセス)
precon_thblock integer in ブロックヤコビ前処理,加法シュワルツ前処理のブロック数
independ_nvec integer in 細分化前処理で独立に計算する対角ブロックのサイズ
nBlock integer in サイクリックにスレッドに割り当てる行数
precision_A integer in 行列の非零要素の精度(1:倍精度,2:4倍精度)
precision_b integer in 右辺ベクトルの精度(1:倍精度,2:4倍精度)
precision_x integer in 反復の初期値,解ベクトルの精度(1:倍精度,2:4倍精度)
precision_precon integer in 前処理行列の精度(1:倍精度,2:4倍精度)
iret integer out エラーコード(0:正常終了)


6.24 parcel_ddgmres_dia


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形式で格納.

引数名(次元) 入力/出力 説明
icomm integer in MPIコミュニケータ
vecx_hi(n) double precision in/out in:反復の初期値(上位ビット),out:解ベクトル(上位ビット)
vecx_lo(n) double precision in/out in:反復の初期値(下位ビット),out:解ベクトル(下位ビット)
vecb_hi(n) double precision in 右辺ベクトル(上位ビット)
vecb_lo(n) double precision in 右辺ベクトル(下位ビット)
n integer*4 / integer*8 in 各プロセスが担当するベクトルの大きさ
gn integer*4 / integer*8 in 全体のベクトルの大きさ
istart integer*4 / integer*8 in 各プロセスが担当する行列の開始行
diaA_hi(n*nnd) double precision in DIA形式に格納された行列の非零要素(上位ビット)
diaA_lo(n*nnd) double precision in DIA形式に格納された行列の非零要素(下位ビット)
offset(nnd) integer*4 / integer*8 in DIA形式に格納された行列の対角からのオフセット
nnd integer*4 / integer*8 in DIA形式に格納された行列の対角要素の本数
ipreflag integer in 前処理指定フラグ(0:無し,1:ポイントヤコビ,2:ブロックヤコビ,3:加法シュワルツ)
ilu_method integer in 不完全LU分解指定フラグ(0:ILU(0),1:D-ILU(対角成分が一致),2:D-ILU(行内の要素和が一致),加法シュワルツは1,2のみ利用可能)
addL integer in 加法シュワルツの重なり幅
iflagAS integer in 加法シュワルツ指定フラグ(1:BASIC,2:RESTRICT,3:INTERPOLATE,4:NONE)
itrmax integer in/out in:最大反復回数,out:反復回数
rtolmax double precision in 収束判定条件(相対残差ノルム)
reshistory(itrmax) double precision out 残差履歴
iovlflag integer in 通信隠蔽フラグ(0:無し,1:全プロセス,2:各プロセス)
gmres_m integer in リスタートするまでの反復回数
gmres_GSflag integer in 直交化フラグ(1:MGS,2:CGS)
precon_thblock integer in ブロックヤコビ前処理,加法シュワルツ前処理のブロック数
independ_nvec integer in 細分化前処理で独立に計算する対角ブロックのサイズ
nBlock integer in サイクリックにスレッドに割り当てる行数
precision_A integer in 行列の非零要素の精度(1:倍精度,2:4倍精度)
precision_b integer in 右辺ベクトルの精度(1:倍精度,2:4倍精度)
precision_x integer in 反復の初期値,解ベクトルの精度(1:倍精度,2:4倍精度)
precision_precon integer in 前処理行列の精度(1:倍精度,2:4倍精度)
iret integer out エラーコード(0:正常終了)


6.25 parcel_ddcg_ddm


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形式で格納.

引数名(次元) 入力/出力 説明
icomm integer in MPIコミュニケータ
x_hi(m) double precision in/out in:反復の初期値(上位ビット),out:解ベクトル(上位ビット)
x_lo(m) double precision in/out in:反復の初期値(下位ビット),out:解ベクトル(下位ビット)
b_hi(m) double precision in 右辺ベクトル(上位ビット)
b_lo(m) double precision in 右辺ベクトル(下位ビット)
m integer in 各プロセスが担当するベクトルの大きさ
nnd integer in プロセスが担当する係数行列要素の非零成分をもつ対角の本数
ioff_dia(nnd) integer in DDM形式に格納された行列の対角からのオフセット
val_dia_hi(m*nnd) double precision in DDM形式に格納された行列の非ゼロ要素(上位ビット)
val_dia_lo(m*nnd) double precision in DDM形式に格納された行列の非ゼロ要素(下位ビット)
num_neighbor_ranks integer in 隣接するプロセス数
neighbor_ranks(num_neighbor_ranks) integer in 隣接するプロセス番号のリスト
margin integer in ioff_gridの絶対値の最大値
ncomp_grids integer in 構造格子の次元数
ndiv_grids(ncomp_grids) integer in 構造格子の各次元の分割数
num_grids(ncomp_grids) integer in プロセスが担当する格子サイズ
ioff_grids(nnd,ncomp_grids) integer in ステンシルの各次元方向の成分
div_direc_th integer in スレッド並列における領域分割の方向(1:x方向,2:y方向,3:z方向)
ipreflag integer in 前処理指定フラグ(0:無し,1:ポイントヤコビ,2:ブロックヤコビ,3:加法シュワルツ
addL integer in 加法シュワルツの重なり幅
iflagAS integer in 加法シュワルツ指定フラグ(1:BASIC,2:RESTRICT,3:INTERPOLATE,4:NONE)
itrmax integer in/out in:最大反復回数,out:収束回数
rtolmax double precision in 収束判定条件(相対残差ノルム)
reshistory(itrmax) double precision out 残差履歴
iovlflag integer in 通信隠蔽フラグ(0:無し,1:全プロセス,2:各プロセス)
precon_thblock integer in ブロックヤコビ前処理,加法シュワルツ前処理のブロック数
independ_nvec integer in 細分化前処理で独立に計算する対角ブロックのサイズ
nBlock integer in サイクリックにスレッドに割り当てる行数
precision_A integer in 行列の非零要素の精度(1:倍精度,2:4倍精度)
precision_b integer in 右辺ベクトルの精度(1:倍精度,2:4倍精度)
precision_x integer in 反復の初期値,解ベクトルの精度(1:倍精度,2:4倍精度)
precision_precon integer in 前処理行列の精度(1:倍精度,2:4倍精度)
iret integer out エラーコード(0:正常終了)


6.26 parcel_ddbicgstab_ddm


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形式で格納.

引数名(次元) 入力/出力 説明
icomm integer in MPIコミュニケータ
x_hi(m) double precision in/out in:反復の初期値(上位ビット),out:解ベクトル(上位ビット)
x_lo(m) double precision in/out in:反復の初期値(下位ビット),out:解ベクトル(下位ビット)
b_hi(m) double precision in 右辺ベクトル(上位ビット)
b_lo(m) double precision in 右辺ベクトル(下位ビット)
m integer in 各プロセスが担当するベクトルの大きさ
nnd integer in プロセスが担当する係数行列要素の非零成分をもつ対角の本数
ioff_dia(nnd) integer in DDM形式に格納された行列の対角からのオフセット
val_dia_hi(m*nnd) double precision in DDM形式に格納された行列の非ゼロ要素(上位ビット)
val_dia_lo(m*nnd) double precision in DDM形式に格納された行列の非ゼロ要素(下位ビット)
num_neighbor_ranks integer in 隣接するプロセス数
neighbor_ranks(num_neighbor_ranks) integer in 隣接するプロセス番号のリスト
margin integer in ioff_gridの絶対値の最大値
ncomp_grids integer in 構造格子の次元数
ndiv_grids(ncomp_grids) integer in 構造格子の各次元の分割数
num_grids(ncomp_grids) integer in プロセスが担当する格子サイズ
ioff_grids(nnd,ncomp_grids) integer in ステンシルの各次元方向の成分
div_direc_th integer in スレッド並列における領域分割の方向(1:x方向,2:y方向,3:z方向)
ipreflag integer in 前処理指定フラグ(0:無し,1:ポイントヤコビ,2:ブロックヤコビ,3:加法シュワルツ
addL integer in 加法シュワルツの重なり幅
iflagAS integer in 加法シュワルツ指定フラグ(1:BASIC,2:RESTRICT,3:INTERPOLATE,4:NONE)
itrmax integer in/out in:最大反復回数,out:収束回数
rtolmax double precision in 収束判定条件(相対残差ノルム)
reshistory(itrmax) double precision out 残差履歴
iovlflag integer in 通信隠蔽フラグ(0:無し,1:全プロセス,2:各プロセス)
precon_thblock integer in ブロックヤコビ前処理,加法シュワルツ前処理のブロック数
independ_nvec integer in 細分化前処理で独立に計算する対角ブロックのサイズ
nBlock integer in サイクリックにスレッドに割り当てる行数
precision_A integer in 行列の非零要素の精度(1:倍精度,2:4倍精度)
precision_b integer in 右辺ベクトルの精度(1:倍精度,2:4倍精度)
precision_x integer in 反復の初期値,解ベクトルの精度(1:倍精度,2:4倍精度)
precision_precon integer in 前処理行列の精度(1:倍精度,2:4倍精度)
iret integer out エラーコード(0:正常終了)


6.27 parcel_ddgmres_ddm


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形式で格納.

引数名(次元) 入力/出力 説明
icomm integer in MPIコミュニケータ
x_hi(m) double precision in/out in:反復の初期値(上位ビット),out:解ベクトル(上位ビット)
x_lo(m) double precision in/out in:反復の初期値(下位ビット),out:解ベクトル(下位ビット)
b_hi(m) double precision in 右辺ベクトル(上位ビット)
b_lo(m) double precision in 右辺ベクトル(下位ビット)
m integer in 各プロセスが担当するベクトルの大きさ
nnd integer in プロセスが担当する係数行列要素の非零成分をもつ対角の本数
ioff_dia(nnd) integer in DDM形式に格納された行列の対角からのオフセット
val_dia_hi(m*nnd) double precision in DDM形式に格納された行列の非ゼロ要素(上位ビット)
val_dia_lo(m*nnd) double precision in DDM形式に格納された行列の非ゼロ要素(下位ビット)
num_neighbor_ranks integer in 隣接するプロセス数
neighbor_ranks(num_neighbor_ranks) integer in 隣接するプロセス番号のリスト
margin integer in ioff_gridの絶対値の最大値
ncomp_grids integer in 構造格子の次元数
ndiv_grids(ncomp_grids) integer in 構造格子の各次元の分割数
num_grids(ncomp_grids) integer in プロセスが担当する格子サイズ
ioff_grids(nnd,ncomp_grids) integer in ステンシルの各次元方向の成分
div_direc_th integer in スレッド並列における領域分割の方向(1:x方向,2:y方向,3:z方向)
ipreflag integer in 前処理指定フラグ(0:無し,1:ポイントヤコビ,2:ブロックヤコビ,3:加法シュワルツ
addL integer in 加法シュワルツの重なり幅
iflagAS integer in 加法シュワルツ指定フラグ(1:BASIC,2:RESTRICT,3:INTERPOLATE,4:NONE)
itrmax integer in/out in:最大反復回数,out:収束回数
rtolmax double precision in 収束判定条件(相対残差ノルム)
reshistory(itrmax) double precision out 残差履歴
iovlflag integer in 通信隠蔽フラグ(0:無し,1:全プロセス,2:各プロセス)
precon_thblock integer in ブロックヤコビ前処理,加法シュワルツ前処理のブロック数
independ_nvec integer in 細分化前処理で独立に計算する対角ブロックのサイズ
nBlock integer in サイクリックにスレッドに割り当てる行数
gmres_m integer in リスタートするまでの反復回数
gmres_GSflag integer in 直交化フラグ(1:MGS,2:CGS)
precision_A integer in 行列の非零要素の精度(1:倍精度,2:4倍精度)
precision_b integer in 右辺ベクトルの精度(1:倍精度,2:4倍精度)
precision_x integer in 反復の初期値,解ベクトルの精度(1:倍精度,2:4倍精度)
precision_precon integer in 前処理行列の精度(1:倍精度,2:4倍精度)
iret integer out エラーコード(0:正常終了)


6.28 parcel_dqr


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


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

引数名(次元) 入力/出力 説明
n integer*4 / integer*8 in 各プロセスが担当する行列の行数
s integer*4 / integer*8 in 行列の列数
icomm integer in MPIコミュニケータ
V(s*n) double precision in QR分解する入力行列 (Column Major)
Q(s*n) double precision out 直交行列 (Column Major)
R(s*s) double precision out 上三角行列 (Column Major)
iQRflag integer in QR分解フラグ(1:MGS,2:CGS,3:TSQR,4:CholeskyQR,5:CholeskyQR2)


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

 PARCELのサブルーチン使用方法についてCG法のFORTRANサンプルコードを例にして説明する.

7.1 CRS形式 (Fortran)

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


  program main
    use mpi
    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 ierr
 
    call MPI_Init(ierr)
 
    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
    call parcel_dcg( &
         MPI_COMM_WORLD, &
         vecx,vecb, &
         n,gn,nnz,istart, &
         crsA,crsRow_ptr,crsCol, &
         ipreflag,ILU_method,addL,iflagAS, &
         itrmax,rtolmax, &
         reshistory, &
         iovlflag,iret &
         precon_thblock, independ_nvec, &
         nblock, &
         iret &
         )
    call MPI_Finalize(ierr)
 
    deallocate(vecx)
    deallocate(vecb)
 
    deallocate(reshistory)
 
  end program main


7.2 DIA形式 (Fortran)

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


  program main
    use mpi
    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 ierr
 
    call MPI_Init(ierr)
 
    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
 
    call parcel_dcg_dia(&
         MPI_COMM_WORLD, &
         vecx,vecb,&
         n,gn,istart,&
         diaA,offset,nnd,&
         ipreflag,ILU_method,addL,iflagAS,&
         itrmax,rtolmax,&
         reshistory,&
         iovlflag,
         precon_thblock, independ_nvec, &
         nblock, &
         iret&
         )
 
    call MPI_Finalize(ierr)
 
    deallocate(reshistory)
 
  end program main


7.3 DDM形式 (Fortran)

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


    program main
      use mpi
      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(:,:)
 
      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)
 
      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 )
 
      call parcel_dcg_ddm(&
           MPI_COMM_WORLD, &
           vecx,vecb,  &
           m, &
           nnd, offset, diaA, &
           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 &
           )
 
      call MPI_Barrier(MPI_COMM_WORLD,ierr)
 
 
      deallocate( diaA )
      deallocate( offset )
      deallocate( offset_dim )
      deallocate( reshistory_DDM )
      deallocate( vecx )
      deallocate( vecx0 )
      deallocate( vecb )
 
      call MPI_Finalize(ierr)
 
    end program main


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

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

8.1 CRS形式 (C言語)

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


 
   #include < stdio.h >
   #include < stdlib.h >
   #include "mpi.h"
   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);
 
     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;
 
     parcel_dcg_(
          &icomm_fort ,
          vecx,
          vecb,
          &n,
          &gn,
          &nnz,
          &istart,
          crsa ,
          crsrow_ptr,
          crscol,
          &ipreflag ,
          &ilu_method ,
          &addl ,
          &iflagas ,
          &max_niter ,
          &rtolmax ,
          reshistory,
          &iovlflag ,
          &precon_thblock ,
          &independ_nvec ,
          &nBlock,
          &iret
          );
     free(vecx);
     free(vecb);
     free(reshistory);
     MPI_Finalize();
     }


8.2 DIA形式 (C言語)

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


 
   #include < stdio.h >
   #include < stdlib.h >
   #include "mpi.h"
   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);
 
     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;
 
     parcel_dcg_dia_(
          &icomm_fort ,
          vecx,
          vecb,
          &n,
          &gn,
          &istart,
          diaA ,
          offset,
          &nnd,
          &ipreflag ,
          &ilu_method ,
          &addl ,
          &iflagas ,
          &max_niter ,
          &rtolmax ,
          reshistory,
          &iovlflag ,
          &precon_thblock ,
          &independ_nvec ,
          &nBlock,
          &iret
          );
     free(vecx);
     free(vecb);
     free(reshistory);
     MPI_Finalize();
 


8.3 DDM形式 (C言語)

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


 
   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);
 
     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);
     parcel_dcg_ddm_(
              &icomm_fort,
              vecx,
              vecb,
              &m,
              &nnd,
              offset,
              diaA,
              &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
              );
 
     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);
     MPI_Finalize();
 
     return 0;
   }



9 参考文献

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