1 クリロフ部分空間法

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

1.1 前処理付共役勾配法

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

  1. 𝐂𝐨𝐦𝐩𝐮𝐭𝐞𝐫0=𝐛𝐀𝐱0𝐟𝐨𝐫𝐬𝐨𝐦𝐞𝐢𝐧𝐢𝐭𝐢𝐚𝐥𝐠𝐮𝐞𝐬𝐬𝐱1=𝐱0\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. 𝐩1=0\mathbf{p}^{\mathbf{- 1}}\mathbf{= 0}
  3. 𝛂1=0\mathbf{\alpha}_{\mathbf{- 1}}\mathbf{= 0\ }\mathbf{\ }
  4. 𝛃1=0\mathbf{\beta}_{\mathbf{- 1}}\mathbf{= 0}
  5. 𝐬𝐨𝐥𝐯𝐞𝐬𝐟𝐫𝐨𝐦𝐊𝐬=𝐫0\mathbf{{solve\ s\ from\ }}\mathbf{K}\mathbf{s = \ }\mathbf{r}^{\mathbf{0}}
  6. 𝛒0=𝐬,𝐫0\mathbf{\rho}_{\mathbf{0}}\mathbf{=}\left\langle \mathbf{s,}\mathbf{r}^{\mathbf{0}} \right\rangle
  7. 𝐟𝐨𝐫𝐢=0,1,2,3\mathbf{for\ i = 0,1,2,3}\mathbf{\ldots}
  8.   𝐩𝐢=𝐬+𝛃𝐢1𝐩𝐢1\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.   𝐱𝐢=𝐱𝐢1+𝛂𝐢1𝐩𝐢1\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.   𝐫𝐢+1=𝐫𝐢𝛂𝐢𝐪𝐢\mathbf{r}^{\mathbf{i + 1}}\mathbf{=}\mathbf{r}^{\mathbf{i}}\mathbf{-}\mathbf{\alpha}_{\mathbf{i}}\mathbf{q}^{\mathbf{i}}
  14.   𝐬𝐨𝐥𝐯𝐞𝐬𝐟𝐫𝐨𝐦𝐊𝐬=𝐫𝐢+1\mathbf{{solve\ s\ from\ }}\mathbf{K}\mathbf{s = \ }\mathbf{r}^{\mathbf{i + 1}}
  15.   𝛒𝐢+1=𝐬,𝐫𝐢+1\mathbf{\rho}_{\mathbf{i + 1}}\mathbf{=}\left\langle \mathbf{s,}\mathbf{r}^{\mathbf{i + 1}} \right\rangle
  16.   𝐢𝐟𝐫𝐢+1/𝐛𝐬𝐦𝐚𝐥𝐥𝐞𝐧𝐨𝐮𝐠𝐡𝐭𝐡𝐞𝐧\mathbf{{if\ }}\left\| \mathbf{r}^{\mathbf{i + 1}} \right\|\mathbf{\ /\ }\left\| \mathbf{b} \right\|\mathbf{{\ small\ enough\ the}}\mathbf{n}
  17.     𝐱𝐢=𝐱𝐢1+𝛂𝐢1𝐩𝐢1\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.   𝛃𝐢=𝛒𝐢+1/𝛒𝐢\mathbf{\beta}_{\mathbf{i}}\mathbf{=}\mathbf{\rho}_{\mathbf{i + 1}}\mathbf{\ /\ }\mathbf{\rho}_{\mathbf{i}}
  21. 𝐞𝐧𝐝\mathbf{{\ en}}\mathbf{d}

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

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

  1. 𝐂𝐨𝐦𝐩𝐮𝐭𝐞𝐫0=𝐛𝐀𝐱0𝐟𝐨𝐫𝐬𝐨𝐦𝐞𝐢𝐧𝐢𝐭𝐢𝐚𝐥𝐠𝐮𝐞𝐬𝐬𝐱1=𝐱0\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. 𝐩0=𝐫0\mathbf{p}_{\mathbf{0}}\mathbf{=}\mathbf{r}_{\mathbf{0}}
  3. 𝐜1=𝐫0,𝐫0\mathbf{c}_{\mathbf{1}}\mathbf{=}\left\langle \mathbf{r}_{\mathbf{0}}\mathbf{,}\mathbf{r}_{\mathbf{0}} \right\rangle
  4. 𝐟𝐨𝐫𝐢=0,1,2,3\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.   𝐜2=𝐫0,𝐪\mathbf{c}_{\mathbf{2}}\mathbf{=}\left\langle \mathbf{r}_{\mathbf{0}}\mathbf{,q} \right\rangle
  8.   𝛂=𝐜1/𝐜2\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.   𝐜3=𝐞,𝐯/𝐯,𝐯\mathbf{c}_{\mathbf{3}}\mathbf{=}\left\langle \mathbf{e,v} \right\rangle\mathbf{\ /\ }\left\langle \mathbf{v,v} \right\rangle
  13.   𝐱𝐢+1=𝐱𝐢+𝛂𝐩̂+𝐜3𝐞̂\mathbf{x}_{\mathbf{i + 1}}\mathbf{=}\mathbf{x}_{\mathbf{i}}\mathbf{+ \alpha}\widehat{\mathbf{p}}\mathbf{+}\mathbf{c}_{\mathbf{3}}\widehat{\mathbf{e}}
  14.   𝐫𝐢+1=𝐞𝐜3𝐯\mathbf{r}_{\mathbf{i + 1}}\mathbf{= e -}\mathbf{c}_{\mathbf{3}}\mathbf{v}
  15.   𝐜1=𝐫0,𝐫𝐢+1\mathbf{c}_{\mathbf{1}}\mathbf{=}\left\langle \mathbf{r}_{\mathbf{0}}\mathbf{,}\mathbf{r}_{\mathbf{i + 1}} \right\rangle
  16.   𝐢𝐟𝐫𝐢+1/𝐛𝐬𝐦𝐚𝐥𝐥𝐞𝐧𝐨𝐮𝐠𝐡𝐭𝐡𝐞𝐧\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.   𝛃=𝐜1/(𝐜2𝐜3)\mathbf{\beta =}\mathbf{c}_{\mathbf{1}}\mathbf{\ /\ }\left( \mathbf{c}_{\mathbf{2}}\mathbf{c}_{\mathbf{3}} \right)
  20.   𝐩𝐢+1=𝐫𝐢+1+𝛃(𝐩𝐢𝐜3𝐪)\mathbf{p}_{\mathbf{i + 1}}\mathbf{=}\mathbf{r}_{\mathbf{i + 1}}\mathbf{+ \beta}\left( \mathbf{p}_{\mathbf{i}}\mathbf{-}\mathbf{c}_{\mathbf{3}}\mathbf{q} \right)
  21. 𝐞𝐧𝐝\mathbf{{end}}

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

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

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

  1. 𝐟𝐨𝐫𝐣=0,1,2,3\mathbf{for\ j = 0,1,2,3}\mathbf{\ldots}
  2.   𝐫=𝐛𝐀𝐱0\mathbf{r = b - A}\mathbf{x}_{\mathbf{0}}
  3.   𝐯1=𝐫/𝐫\mathbf{v}_{\mathbf{1}}\mathbf{= - r\ /\ }\left\| \mathbf{r} \right\|
  4.   𝐞=(𝐫,0,,0)𝐓\mathbf{e =}\left( \mathbf{-}\left\| \mathbf{r} \right\|\mathbf{,0,\ldots,0} \right)^{\mathbf{T}}
  5.   𝐧=𝐦\mathbf{n = m}
  6.   𝐟𝐨𝐫𝐢=1,2,3𝐦\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.      𝐟𝐨𝐫𝐤=1,2,3𝐢\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.     𝐡𝐢+1,𝐢=𝛚\mathbf{h}_{\mathbf{i + 1,i}}\mathbf{=}\left\| \mathbf{\omega} \right\|
  14.     𝐯𝐢+1=𝛚/𝛚\mathbf{v}_{\mathbf{i + 1}}\mathbf{= \omega/}\left\| \mathbf{\omega} \right\|
  15.     𝐟𝐨𝐫𝐤=1,2,3𝐢1\mathbf{for\ k = 1,2,3}\mathbf{\ldots}\mathbf{i - 1}
  16.       (𝐡𝐤,𝐢𝐡𝐤+1,𝐢)=(𝐜𝐢𝐬𝐢𝐬𝐢𝐜𝐢)(𝐡𝐤,𝐢𝐡𝐤+1,𝐢)\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.     𝐜𝐢=11+(𝐡𝐢+1,𝐢𝐡𝐢,𝐢)2\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.     𝐬𝐢=𝐡𝐢+1,𝐢𝐡𝐢,𝐢11+(𝐡𝐢+1,𝐢𝐡𝐢,𝐢)2\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.     (𝐞𝐢𝐞𝐢+1)=(𝐜𝐢𝐬𝐢𝐬𝐢𝐜𝐢)(𝐞𝐢𝐞𝐢+1)\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.     𝐡𝐢,𝐢=𝐜𝐢𝐡𝐢,𝐢𝐬𝐢𝐡𝐢+1,𝐢\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.     𝐡𝐢+1,𝐢=0\mathbf{h}_{\mathbf{i + 1,i}}\mathbf{= 0}
  23.     𝐢𝐟𝐞𝐢+1/𝐛𝐬𝐦𝐚𝐥𝐥𝐞𝐧𝐨𝐮𝐠𝐡𝐭𝐡𝐞𝐧\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.   𝐲𝐧=𝐇𝐧1𝐞𝐧\mathbf{y}_{\mathbf{n}}\mathbf{=}\mathbf{H}_{\mathbf{n}}^{\mathbf{- 1}}\mathbf{e}_{\mathbf{n}}
  29.   𝐬𝐨𝐥𝐯𝐞𝐱̂𝐟𝐫𝐨𝐦𝐊𝐱̂=𝐤=1𝐧𝐲𝐤𝐯𝐤\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.   𝐱𝐧=𝐱0+𝐱̂\mathbf{x}_{\mathbf{n}}\mathbf{=}\mathbf{x}_{\mathbf{0}}\mathbf{+}\widehat{\mathbf{x}}
  31.   𝐱0=𝐱𝐧\mathbf{x}_{\mathbf{0}}\mathbf{=}\mathbf{x}_{\mathbf{n}}
  32. 𝐞𝐧𝐝\mathbf{{end}}

1.4 前処理付通信削減型一般化最小残差法

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

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

  1. 𝐁𝐤=[𝐞2,𝐞3,,𝐞𝐬+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. 𝐟𝐨𝐫𝐣=0,1,2,3\mathbf{for\ j = 0,1,2,3\ldots}
  3.   𝐫=𝐛𝐀𝐱0\mathbf{r = b - A}\mathbf{x}_{\mathbf{0}}
  4.   𝐯1=𝐫/𝐫\mathbf{v}_{\mathbf{1}}\mathbf{= r\ /\ }\left\| \mathbf{r} \right\|
  5.   𝛇=(𝐫,0,,0)𝐓\mathbf{\zeta =}\left( \left\| \mathbf{r} \right\|\mathbf{,0,\ldots,0} \right)^{\mathbf{T}}
  6.   𝐟𝐨𝐫𝐤=0,1,𝐭1\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.     𝐢𝐟(𝐤.𝐞𝐪.0)\mathbf{if(k.eq.0)}
  10.       𝐐𝐑𝐝𝐞𝐜𝐨𝐦𝐩𝐨𝐬𝐢𝐭𝐢𝐨𝐧𝐕0=𝐐0𝐑0\mathbf{{QR\ decomposition\ }}\mathbf{V}_{\mathbf{0}}^{\mathbf{\sim}}\mathbf{=}\mathbf{Q}_{\mathbf{0}}^{\mathbf{\sim}}\mathbf{R}_{\mathbf{0}}^{\mathbf{\sim}}
  11.       𝔔0=𝐐0\mathfrak{Q}_{\mathbf{0}}^{\mathbf{\sim}}\mathbf{=}\mathbf{Q}_{\mathbf{0}}^{\mathbf{\sim}}
  12.       0=𝐑0𝐁0𝐑01\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.       0=0\mathcal{H}_{\mathbf{0}}\mathbf{=}\mathfrak{H}_{\mathbf{0}}^{\mathbf{\sim}}
  14.       𝐟𝐨𝐫𝐨=1,2,𝐬\mathbf{for\ o = 1,2\ldots,s}
  15.       𝐆𝐢𝐯𝐞𝐧𝐬𝐫𝐨𝐭𝐚𝐭𝐢𝐨𝐧[𝐨,0,𝛇]\mathbf{Givens\ rotation\ \lbrack o,}\mathcal{H}_{\mathbf{0}}\mathbf{,\zeta\rbrack}
  16.       𝐢𝐟𝛇𝐨+1/𝐛𝐬𝐦𝐚𝐥𝐥𝐞𝐧𝐨𝐮𝐠𝐡𝐭𝐡𝐞𝐧\mathbf{{if\ }}\left\| \mathbf{\zeta}_{\mathbf{o + 1}} \right\|\mathbf{\ /\ }\left\| \mathbf{b} \right\|\mathbf{{\ small\ enough\ then}}
  17.         𝐮𝐩𝐝𝐚𝐭𝐞𝐬𝐨𝐥𝐮𝐭𝐢𝐨𝐧𝐯𝐞𝐜𝐭𝐨𝐫[𝐬,𝐤,𝐨,𝐕́0,𝐐́0,𝐑0,0,𝛇,𝐱0]\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.       ́𝐤1,𝐤=(𝔔𝐤1)𝐇𝐕́𝐤{\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.       𝐕́𝐤=𝐕́𝐤𝔔𝐤1́𝐤1,𝐤{\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.       𝐑𝐤=(𝐑𝐤𝐳𝐤01,𝐤𝛒)\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.       𝐤1,𝐤=𝐤1𝐤1,𝐤𝐑𝐤1+𝐤1,𝐤𝐁𝐤𝐑𝐤1\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.       𝐇𝐤=𝐑𝐤𝐁𝐤𝐑𝐤1+𝛒̃𝐤1𝐛𝐤𝐳𝐤𝐞𝐬𝐓𝐡𝐤1𝐞1𝐞𝐬(𝐤1)𝐓𝐤1,𝐤𝐑𝐤1\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.       𝐡𝐤=𝛒̃𝐤1𝛒𝐤𝐛𝐤\mathbf{h}_{\mathbf{k}}\mathbf{=}{\widetilde{\mathbf{\rho}}}_{\mathbf{k}}^{\mathbf{- 1}}\mathbf{\rho}_{\mathbf{k}}\mathbf{b}_{\mathbf{k}}
  28.       𝐤=(𝐤1𝐤1,𝐤𝐡𝐤1𝐞1𝐞𝐬(𝐤1)𝐓01,𝐬(𝐤1)𝐇𝐤𝐡𝐤𝐞𝐬𝐓)\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.       𝐤=(𝐤1,𝐤𝐇𝐤𝐡𝐤𝐞𝐬𝐓)\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.       𝐟𝐨𝐫𝐨=1+𝐬𝐤,,𝐬(𝐤+1)\mathbf{for\ o = 1 + sk,\ldots,s(k + 1)}
  31.         𝐆𝐢𝐯𝐞𝐧𝐬𝐫𝐨𝐭𝐚𝐭𝐢𝐨𝐧[𝐨,𝐤,𝛇]\mathbf{Givens\ rotation\ \lbrack o,}\mathcal{H}_{\mathbf{k}}\mathbf{,\zeta\rbrack}
  32.         𝐢𝐟𝛇𝐨+1/𝐛𝐬𝐦𝐚𝐥𝐥𝐞𝐧𝐨𝐮𝐠𝐡𝐭𝐡𝐞𝐧\mathbf{{if\ }}\left\| \mathbf{\zeta}_{\mathbf{o + 1}} \right\|\mathbf{\ /\ }\left\| \mathbf{b} \right\|\mathbf{{\ small\ enough\ then}}
  33.           𝐮𝐩𝐝𝐚𝐭𝐞𝐬𝐨𝐥𝐮𝐭𝐢𝐨𝐧𝐯𝐞𝐜𝐭𝐨𝐫[𝐬,𝐤,𝐨,𝐕́𝐤,𝔔𝐤,𝐑́𝐤,𝐤,𝛇,𝐱0]\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.   𝐮𝐩𝐝𝐚𝐭𝐞𝐬𝐨𝐥𝐮𝐭𝐢𝐨𝐧𝐯𝐞𝐜𝐭𝐨𝐫[𝐬,𝐭1,𝐬𝐭,𝐕́𝐭1,𝔔𝐭1,𝐑́𝐭1,𝐭1,𝛇,𝐱0]\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.   𝐢=0\mathbf{i = 0}
  3.   𝐰𝐡𝐢𝐥𝐞(𝐢𝐬1)\mathbf{{while}}\left( \mathbf{i \leq s - 1} \right)
  4.     𝐢𝐟(𝐢.𝐞𝐪.𝐬1)𝐭𝐡𝐞𝐧\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.       𝐢𝐟(𝐈𝐦[𝐄𝐢].𝐞𝐪.0)𝐭𝐡𝐞𝐧\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.         𝐁𝐢+1,𝐢+1=𝐑𝐞[𝐄𝐢]\mathbf{B}_{\mathbf{i + 1,i + 1}}\mathbf{= Re\lbrack}\mathbf{E}_{\mathbf{i}}\mathbf{\rbrack}
  12.         𝐁𝐢,𝐢+1=(𝐈𝐦[𝐄𝐢])2\mathbf{B}_{\mathbf{i,i + 1}}\mathbf{= -}\left( \mathbf{Im\lbrack}\mathbf{E}_{\mathbf{i}}\mathbf{\rbrack} \right)^{\mathbf{2}}
  13.         𝐢=𝐢+1\mathbf{i = i + 1}
  14.       𝐞𝐧𝐝𝐢𝐟\mathbf{{endif}}
  15.     𝐞𝐧𝐝𝐢𝐟\mathbf{{endif}}
  16.     𝐢=𝐢+1\mathbf{i = i + 1}
  17.   𝐞𝐧𝐝\mathbf{{end\ \ }}
  18. 𝐞𝐧𝐝\mathbf{{end\ \ }}

𝐁\mathbf{B}𝐛𝐤,𝐢\mathbf{b}_{\mathbf{k,i}}を要素に持つ行列.

  1. 𝐟𝐨𝐫𝐤=1,2,3𝐬\mathbf{for\ k = 1,2,3\ldots s}
  2.   𝐬𝐨𝐥𝐯𝐞𝐯̂𝐤1𝐟𝐫𝐨𝐦𝐊𝐯̂𝐤1=𝐯𝐤1\mathbf{{solve\ }}{\widehat{\mathbf{v}}}_{\mathbf{k - 1}}\mathbf{{\ from\ K}}{\widehat{\mathbf{v}}}_{\mathbf{k - 1}}\mathbf{=}\mathbf{v}_{\mathbf{k - 1}}
  3.   𝐢𝐟(𝐤.𝐧𝐞.1)𝐭𝐡𝐞𝐧\mathbf{if(k.ne.1)then}
  4.     𝛂=𝐛𝐤1,𝐤1\mathbf{\alpha =}\mathbf{b}_{\mathbf{k - 1,k - 1}}
  5.     𝛃=𝐛𝐤2,𝐤1\mathbf{\beta =}\mathbf{b}_{\mathbf{k - 2,k - 1}}
  6.     𝐯𝐤=𝐀𝐯̂𝐤1𝛂𝐯𝐤1+𝛃𝐯𝐤2\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.     𝛂=𝐛𝐤1,𝐤1\mathbf{\alpha =}\mathbf{b}_{\mathbf{k - 1,k - 1}}
  9.     𝐯𝐤=𝐀𝐯̂𝐤1𝛂𝐯𝐤1\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. 𝐟𝐨𝐫𝐤=1,2,3𝐢1\mathbf{for\ k = 1,2,3\ldots i - 1}
  2.   (𝒽𝐤,𝐢𝒽𝐤+1,𝐢)=(𝐜𝐢𝐬𝐢𝐬𝐢𝐜𝐢)(𝒽𝐤,𝐢𝒽𝐤+1,𝐢)\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. end\mathbf{\text{end}}
  4. 𝐜𝐢=11+(𝒽𝐢+1,𝐢𝒽𝐢,𝐢)2\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. 𝐬𝐢=𝒽𝐢+1,𝐢𝒽𝐢,𝐢11+(𝒽𝐢+1,𝐢𝒽𝐢,𝐢)2\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. (𝛇𝐢𝛇𝐢+1)=(𝐜𝐢𝐬𝐢𝐬𝐢𝐜𝐢)(𝛇𝐢𝛇𝐢+1)\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. 𝒽𝐢,𝐢=𝐜𝐢𝒽𝐢,𝐢𝐬𝐢𝒽𝐢+1,𝐢\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. 𝒽𝐢+1,𝐢=0\mathcal{h}_{\mathbf{i + 1,i}}\mathbf{= 0}

  1. 𝐲=1𝛇\mathbf{y =}\mathcal{H}^{\mathbf{- 1}}\mathbf{\zeta}
  2. 𝐬𝐨𝐥𝐯𝐞𝐱̂𝐟𝐫𝐨𝐦𝐊𝐱̂=𝐤=0𝐬𝐦1𝐐𝐤𝐲𝐤+𝐥=𝐬𝐦𝐧1𝐕𝐥𝐬𝐦𝐤=𝐬𝐦𝐧1𝐑𝐥𝐬𝐦,𝐤𝐬𝐦1𝐲𝐤\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. 𝐱0=𝐱0+𝐱̂\mathbf{x}_{\mathbf{0}}\mathbf{=}\mathbf{x}_{\mathbf{0}}\mathbf{+}\widehat{\mathbf{x}}

1.5 前処理付チェビシェフ基底共役勾配法法

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

  1. 𝐂𝐨𝐦𝐩𝐮𝐭𝐞𝐫0=𝐛𝐀𝐱0𝐟𝐨𝐫𝐬𝐨𝐦𝐞𝐢𝐧𝐢𝐭𝐢𝐚𝐥𝐠𝐮𝐞𝐬𝐬𝐱1=𝐱0\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. 𝐂𝐨𝐦𝐩𝐮𝐭𝐞𝐂𝐡𝐞𝐛𝐲𝐬𝐡𝐞𝐯𝐛𝐚𝐬𝐢𝐬[𝐒0,𝐀𝐒0,𝐫0,𝛌max,𝛌min]\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. 𝐐0=𝐒0\mathbf{Q}_{\mathbf{0}}\mathbf{=}\mathbf{S}_{\mathbf{0}}
  4. 𝐟𝐨𝐫𝐢=0,1,2,3\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.   𝛂𝐢=(𝐐𝐢𝐓𝐀𝐐𝐢)1(𝐐𝐢𝐓𝐫𝐢𝐤)\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.   𝐱(𝐢+1)𝐤=𝐱𝐢𝐤+𝐐𝐢𝛂𝐢\mathbf{x}_{\left( \mathbf{i + 1} \right)\mathbf{k}}\mathbf{=}\mathbf{x}_{\mathbf{{ik}}}\mathbf{+}\mathbf{Q}_{\mathbf{i}}\mathbf{\alpha}_{\mathbf{i}}
  8.   𝐫(𝐢+1)𝐤=𝐫𝐢𝐤𝐀𝐐𝐢𝛂𝐢\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.   𝐢𝐟𝐫(𝐢+1)𝐤/𝐛𝐬𝐦𝐚𝐥𝐥𝐞𝐧𝐨𝐮𝐠𝐡𝐭𝐡𝐞𝐧\mathbf{{if\ }}\left\| \mathbf{r}_{\left( \mathbf{i + 1} \right)\mathbf{k}} \right\|\mathbf{\ /\ }\left\| \mathbf{b} \right\|\mathbf{{\ small\ enough\ }}\mathbf{{then}}
  10.   𝐂𝐨𝐦𝐩𝐮𝐭𝐞𝐂𝐡𝐞𝐛𝐲𝐬𝐡𝐞𝐯𝐛𝐚𝐬𝐢𝐬[𝐒𝐢+1,𝐀𝐒𝐢+1,𝐫(𝐢+1)𝐤,𝛌max,𝛌min]\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.   𝐂𝐨𝐦𝐩𝐮𝐭𝐞𝐐𝐢𝐓𝐀𝐒𝐢+1\mathbf{{Compute\ }}\mathbf{Q}_{\mathbf{i}}^{\mathbf{T}}\mathbf{A}\mathbf{S}_{\mathbf{i + 1}}
  12.   𝐁𝐢=(𝐐𝐢𝐓𝐀𝐐𝐢)1(𝐐𝐢𝐓𝐀𝐒𝐢+1)\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.   𝐐𝐢+1=𝐒𝐢+1𝐐𝐢𝐁𝐢\mathbf{Q}_{\mathbf{i + 1}}\mathbf{=}\mathbf{S}_{\mathbf{i + 1}}\mathbf{-}\mathbf{Q}_{\mathbf{i}}\mathbf{B}_{\mathbf{i}}
  14.   𝐀𝐐𝐢+1=𝐀𝐒𝐢+1𝐀𝐐𝐢𝐁𝐢\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}}

𝛌max\mathbf{\lambda}_{\mathbf{\max}},𝛌min\mathbf{\lambda}_{\mathbf{\min}}𝐀𝐊1\mathbf{AK}^{-1}の最大及び最小固有値.

  1. 𝛈=2/(𝛌max𝛌min)\mathbf{\eta = 2\ /\ (}\mathbf{\lambda}_{\mathbf{\max}}\mathbf{-}\mathbf{\lambda}_{\mathbf{\min}}\mathbf{)}
  2. 𝛇=(𝛌max+𝛌min)/(𝛌max𝛌min)\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. 𝐬0=𝐫0\mathbf{s}_{\mathbf{0}}\mathbf{=}\mathbf{r}_{\mathbf{0}}\mathbf{\ }
  4. 𝐬𝐨𝐥𝐯𝐞𝐬̃0𝐟𝐫𝐨𝐦𝐊𝐬̃0=𝐬0\mathbf{{solve\ }}{\widetilde{\mathbf{s}}}_{\mathbf{0}}\mathbf{{\ from\ K}}{\widetilde{\mathbf{s}}}_{\mathbf{0}}\mathbf{= \ }\mathbf{s}_{\mathbf{0}}
  5. 𝐬1=𝛈𝐀𝐬̃0𝛇𝐬0\mathbf{s}_{\mathbf{1}}\mathbf{= \eta A}{\widetilde{\mathbf{s}}}_{\mathbf{0}}\mathbf{- \zeta}\mathbf{s}_{\mathbf{0}}
  6. 𝐬𝐨𝐥𝐯𝐞𝐬̃1𝐟𝐫𝐨𝐦𝐊𝐬̃1=𝐬1\mathbf{{solve\ }}{\widetilde{\mathbf{s}}}_{\mathbf{1}}\mathbf{{\ from\ K}}{\widetilde{\mathbf{s}}}_{\mathbf{1}}\mathbf{= \ }\mathbf{s}_{\mathbf{1}}
  7. 𝐟𝐨𝐫𝐣=2,3,,𝐤\mathbf{for\ j = 2,3,\ldots,k}
  8.   𝐬𝐣=2𝛈𝐀𝐬̃𝐣2𝛇𝐬𝐣1𝐬𝐣2\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. 𝐒=(𝐬̃0,𝐬̃1,,𝐬̃𝐤1)\mathbf{S = (}{\widetilde{\mathbf{s}}}_{\mathbf{0}}\mathbf{,}{\widetilde{\mathbf{s}}}_{\mathbf{1}}\mathbf{,\ldots,}{\widetilde{\mathbf{s}}}_{\mathbf{k - 1}}\mathbf{)}
  12. 𝐀𝐒=(𝐀𝐬̃0,𝐀𝐬̃1,,𝐀𝐬̃𝐤1)\mathbf{AS = (A}{\widetilde{\mathbf{s}}}_{\mathbf{0}}\mathbf{,}{\mathbf{A}\widetilde{\mathbf{s}}}_{\mathbf{1}}\mathbf{,\ldots,}{\mathbf{A}\widetilde{\mathbf{s}}}_{\mathbf{k - 1}}\mathbf{)}

2 前処理

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

2.1 ポイントヤコビ前処理

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

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

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

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

 係数行列AAをその対角行列DD,下三角行列LL,上三角行列UUを用いて A=LA+DA+UAA = L_{A} + D_{A} + U_{A} と分解する.このLAL_{A},UAU_{A}および対角行列DD (DDAD \neq D_{A})を用いて,前処理行列KKK=(D+LA)D1(D+UA)K = \left( D + L_{A} \right)D^{- 1}\left( D + U_{A} \right) の形に決定する.ここでは,対角行列DDの決定方法として,次の2とおりの条件のどちらを選択する.

対角行列DDの具体的な計算方法は

となる.

2.4 ブロックヤコビ前処理

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

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

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

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

BASIC

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

図2 BASIC

RESTRICT

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

図3 RESTRICT

INTERPOLATE

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

図4 INTERPOLATE

NONE

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

図5 NONE

3 データ構造

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

3.1 Compressed Row 形式

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


A=(abc00d00e0fg0hij) 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}crsA=\{a,b,c,d,e,f,g,h,i,j\}

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

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


2プロセスの場合

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

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

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

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

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

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


3.2 Diagonal形式

 Diagonal 形式(圧縮対角形式)は,対角行単位に格納する方法である.対角行の値と対角からのオフセットを1次元配列に格納する. 並列処理の場合には列方向に行列を分割して格納する.


A=(abc00d00e0fg0hij) 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}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}offset=\{-2,-1,0,1,2\}

対角要素の本数 :: nnd=5nnd=5


2プロセスの場合

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

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

対角要素の本数 :: nnd=3nnd=3

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

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

対角要素の本数 :: nnd=4nnd=4


4 PARCELの構成

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

4.1 ディレクトリ構成

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

図6 PARCELディレクトリ構成

4.2 コンパイル方法

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


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

図7 make_configの記述例

4.3 使用方法

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

5 PARCELサブルーチン(直接インターフェース)

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

5.1 parcel_dcg


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


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

引数名(次元) 入力/出力 説明
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 Compressed Row 形式に格納された行列の非零要素
crsRow_ptr(n+1) integer*4 / integer*8 in Compressed Row 形式のポインタテーブル
crsCol(nnz) integer*4 / integer*8 in Compressed Row 形式に格納された分割行列の非零要素の列番号
ipreflag integer in 前処理指定フラグ(0:無し,1:ポイントヤコビ,2:ブロックヤコビ,3:加法シュワルツ)
ILU_method integer in 不完全LU分解指定フラグ(0:ILU(0),1:D-ILU(対角成分が一致),2:D-ILU(行内の要素和が一致)
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:各プロセス)
iret integer out エラーコード(0:正常終了)

5.2 parcel_dcg_dia


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


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

引数名(次元) 入力/出力 説明
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 Diagonal形式に格納された行列の非零要素
offset(nnd) integer*4 / integer*8 in Diagonal形式に格納された行列の対角からのオフセット
nnd integer*4 / integer*8 in Diagonal形式に格納された行列の対角要素の本数
ipreflag integer in 前処理指定フラグ(0:無し,1:ポイントヤコビ,2:ブロックヤコビ,3:加法シュワルツ)
ILU_method integer in 不完全LU分解指定フラグ(0:ILU(0),1:D-ILU(対角成分が一致),2:D-ILU(行内の要素和が一致)
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:各プロセス)
iret integer out エラーコード(0:正常終了)

5.3 parcel_ddcg


call parcel_ddcg( icomm, vecx_hi,vecx_lo, vecb_hi,vecb_lo, n,gn,nnz,istart, crsA_hi,crsA_lo,crsRow_ptr,crsCol, ipreflag,ILU_method,addL,iflagAS, itrmax,rtolmax, reshistory, iovlflag,iret, PRECISION_A, PRECISION_b, PRECISION_x, PRECISION_PRECONDITION )


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

引数名(次元) 入力/出力 説明
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 Compressed Row 形式に格納された行列の非零要素(上位ビット)
crsA_lo(nnz) double precision in Compressed Row 形式に格納された行列の非零要素(下位ビット)
crsRow_ptr(n+1) integer*4 / integer*8 in Compressed Row 形式のポインタテーブル
crsCol(nnz) integer*4 / integer*8 in Compressed Row 形式に格納された分割行列の非零要素の列番号
ipreflag integer in 前処理指定フラグ(0:無し,1:ポイントヤコビ,2:ブロックヤコビ,3:加法シュワルツ)
ILU_method integer in 不完全LU分解指定フラグ(0:ILU(0),1:D-ILU(対角成分が一致),2:D-ILU(行内の要素和が一致)
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:各プロセス)
iret integer out エラーコード(0:正常終了)
PRECISION_A integer in 行列の非零要素の精度(1:倍精度,2:4倍精度)
PRECISION_b integer in 右辺ベクトルの精度(1:倍精度,2:4倍精度)
PRECISION_x integer in 反復の初期値,解ベクトルの精度(1:倍精度,2:4倍精度)
PRECISION_PRECONDITION integer in 前処理行列の精度(1:倍精度,2:4倍精度)

5.4 parcel_ddcg_dia


call parcel_ddcg_dia( icomm, vecx_hi,vecx_lo, vecb_hi,vecb_lo, n,gn,istart, diaA_hi,diaA_lo,offset,nnd, ipreflag,ILU_method,addL,iflagAS, itrmax,rtolmax, reshistory, iovlflag,iret, PRECISION_A, PRECISION_b, PRECISION_x, PRECISION_PRECONDITION )


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

引数名(次元) 入力/出力 説明
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 Diagonal形式に格納された行列の非零要素(上位ビット)
diaA_lo(n*nnd) double precision in Diagonal形式に格納された行列の非零要素(下位ビット)
offset(nnd) integer*4 / integer*8 in Diagonal形式に格納された行列の対角からのオフセット
nnd integer*4 / integer*8 in Diagonal形式に格納された行列の対角要素の本数
ipreflag integer in 前処理指定フラグ(0:無し,1:ポイントヤコビ,2:ブロックヤコビ,3:加法シュワルツ)
ILU_method integer in 不完全LU分解指定フラグ(0:ILU(0),1:D-ILU(対角成分が一致),2:D-ILU(行内の要素和が一致)
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:各プロセス)
iret integer out エラーコード(0:正常終了)
PRECISION_A integer in 行列の非零要素の精度(1:倍精度,2:4倍精度)
PRECISION_b integer in 右辺ベクトルの精度(1:倍精度,2:4倍精度)
PRECISION_x integer in 反復の初期値,解ベクトルの精度(1:倍精度,2:4倍精度)
PRECISION_PRECONDITION integer in 前処理行列の精度(1:倍精度,2:4倍精度)

5.5 parcel_dbicgstab


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


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

引数名(次元) 入力/出力 説明
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 Compressed Row 形式に格納された行列の非零要素
crsRow_ptr(n+1) integer*4 / integer*8 in Compressed Row 形式のポインタテーブル
crsCol(nnz) integer*4 / integer*8 in Compressed Row 形式に格納された分割行列の非零要素の列番号
ipreflag integer in 前処理指定フラグ(0:無し,1:ポイントヤコビ,2:ブロックヤコビ,3:加法シュワルツ)
ILU_method integer in 不完全LU分解指定フラグ(0:ILU(0),1:D-ILU(対角成分が一致),2:D-ILU(行内の要素和が一致)
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:各プロセス)
iret integer out エラーコード(0:正常終了)

5.6 parcel_dbicgstab_dia


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


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

引数名(次元) 入力/出力 説明
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 Diagonal形式に格納された行列の非零要素
offset(nnd) integer*4 / integer*8 in Diagonal形式に格納された行列の対角からのオフセット
nnd integer*4 / integer*8 in Diagonal形式に格納された行列の対角要素の本数
ipreflag integer in 前処理指定フラグ(0:無し,1:ポイントヤコビ,2:ブロックヤコビ,3:加法シュワルツ)
ILU_method integer in 不完全LU分解指定フラグ(0:ILU(0),1:D-ILU(対角成分が一致),2:D-ILU(行内の要素和が一致)
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:各プロセス)
iret integer out エラーコード(0:正常終了)

5.7 parcel_ddbicgstab


call parcel_ddbicgstab( icomm, vecx_hi,vecx_lo, vecb_hi,vecb_lo, n,gn,nnz,istart, crsA_hi,crsA_lo,crsRow_ptr,crsCol, ipreflag,ILU_method,addL,iflagAS, itrmax,rtolmax, reshistory, iovlflag,iret, PRECISION_A, PRECISION_b, PRECISION_x, PRECISION_PRECONDITION )


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

引数名(次元) 入力/出力 説明
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 Compressed Row 形式に格納された行列の非零要素(上位ビット)
crsA_lo(nnz) double precision in Compressed Row 形式に格納された行列の非零要素(下位ビット)
crsRow_ptr(n+1) integer*4 / integer*8 in Compressed Row 形式のポインタテーブル
crsCol(nnz) integer*4 / integer*8 in Compressed Row 形式に格納された分割行列の非零要素の列番号
ipreflag integer in 前処理指定フラグ(0:無し,1:ポイントヤコビ,2:ブロックヤコビ,3:加法シュワルツ)
ILU_method integer in 不完全LU分解指定フラグ(0:ILU(0),1:D-ILU(対角成分が一致),2:D-ILU(行内の要素和が一致)
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:各プロセス)
iret integer out エラーコード(0:正常終了)
PRECISION_A integer in 行列の非零要素の精度(1:倍精度,2:4倍精度)
PRECISION_b integer in 右辺ベクトルの精度(1:倍精度,2:4倍精度)
PRECISION_x integer in 反復の初期値,解ベクトルの精度(1:倍精度,2:4倍精度)
PRECISION_PRECONDITION integer in 前処理行列の精度(1:倍精度,2:4倍精度)

5.8 parcel_ddbicgstab_dia


call parcel_ddbicgstab_dia( icomm, vecx_hi,vecx_lo, vecb_hi,vecb_lo, n,gn,istart, diaA_hi,diaA_lo,offset,nnd, ipreflag,ILU_method,addL,iflagAS, itrmax,rtolmax, reshistory, iovlflag,iret, PRECISION_A, PRECISION_b, PRECISION_x, PRECISION_PRECONDITION )


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

引数名(次元) 入力/出力 説明
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 Diagonal形式に格納された行列の非零要素(上位ビット)
diaA_lo(n*nnd) double precision in Diagonal形式に格納された行列の非零要素(下位ビット)
offset(nnd) integer*4 / integer*8 in Diagonal形式に格納された行列の対角からのオフセット
nnd integer*4 / integer*8 in Diagonal形式に格納された行列の対角要素の本数
ipreflag integer in 前処理指定フラグ(0:無し,1:ポイントヤコビ,2:ブロックヤコビ,3:加法シュワルツ)
ILU_method integer in 不完全LU分解指定フラグ(0:ILU(0),1:D-ILU(対角成分が一致),2:D-ILU(行内の要素和が一致)
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:各プロセス)
iret integer out エラーコード(0:正常終了)
PRECISION_A integer in 行列の非零要素の精度(1:倍精度,2:4倍精度)
PRECISION_b integer in 右辺ベクトルの精度(1:倍精度,2:4倍精度)
PRECISION_x integer in 反復の初期値,解ベクトルの精度(1:倍精度,2:4倍精度)
PRECISION_PRECONDITION integer in 前処理行列の精度(1:倍精度,2:4倍精度)

5.9 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,iret, GMRES_m,GMRES_GSflag )


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

引数名(次元) 入力/出力 説明
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 Compressed Row 形式に格納された行列の非零要素
crsRow_ptr(n+1) integer*4 / integer*8 in Compressed Row 形式のポインタテーブル
crsCol(nnz) integer*4 / integer*8 in Compressed Row 形式に格納された分割行列の非零要素の列番号
ipreflag integer in 前処理指定フラグ(0:無し,1:ポイントヤコビ,2:ブロックヤコビ,3:加法シュワルツ)
ILU_method integer in 不完全LU分解指定フラグ(0:ILU(0),1:D-ILU(対角成分が一致),2:D-ILU(行内の要素和が一致)
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:各プロセス)
iret integer out エラーコード(0:正常終了)
GMRES_m integer in リスタートするまでの反復回数
GMRES_GSflag integer in 直交化フラグ(1:MGS,2:CGS)

5.10 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,iret, GMRES_m,GMRES_GSflag )


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

引数名(次元) 入力/出力 説明
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 Diagonal形式に格納された行列の非零要素
offset(nnd) integer*4 / integer*8 in Diagonal形式に格納された行列の対角からのオフセット
nnd integer*4 / integer*8 in Diagonal形式に格納された行列の対角要素の本数
ipreflag integer in 前処理指定フラグ(0:無し,1:ポイントヤコビ,2:ブロックヤコビ,3:加法シュワルツ)
ILU_method integer in 不完全LU分解指定フラグ(0:ILU(0),1:D-ILU(対角成分が一致),2:D-ILU(行内の要素和が一致)
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:各プロセス)
iret integer out エラーコード(0:正常終了)
GMRES_m integer in リスタートするまでの反復回数
GMRES_GSflag integer in 直交化フラグ(1:MGS,2:CGS)

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


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

引数名(次元) 入力/出力 説明
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 Compressed Row 形式に格納された行列の非零要素(上位ビット)
crsA_lo(nnz) double precision in Compressed Row 形式に格納された行列の非零要素(下位ビット)
crsRow_ptr(n+1) integer*4 / integer*8 in Compressed Row 形式のポインタテーブル
crsCol(nnz) integer*4 / integer*8 in Compressed Row 形式に格納された分割行列の非零要素の列番号
ipreflag integer in 前処理指定フラグ(0:無し,1:ポイントヤコビ,2:ブロックヤコビ,3:加法シュワルツ)
ILU_method integer in 不完全LU分解指定フラグ(0:ILU(0),1:D-ILU(対角成分が一致),2:D-ILU(行内の要素和が一致)
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:各プロセス)
iret integer out エラーコード(0:正常終了)
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_PRECONDITION integer in 前処理行列の精度(1:倍精度,2:4倍精度)

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


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

引数名(次元) 入力/出力 説明
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 Diagonal形式に格納された行列の非零要素(上位ビット)
diaA_lo(n*nnd) double precision in Diagonal形式に格納された行列の非零要素(下位ビット)
offset(nnd) integer*4 / integer*8 in Diagonal形式に格納された行列の対角からのオフセット
nnd integer*4 / integer*8 in Diagonal形式に格納された行列の対角要素の本数
ipreflag integer in 前処理指定フラグ(0:無し,1:ポイントヤコビ,2:ブロックヤコビ,3:加法シュワルツ)
ILU_method integer in 不完全LU分解指定フラグ(0:ILU(0),1:D-ILU(対角成分が一致),2:D-ILU(行内の要素和が一致)
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:各プロセス)
iret integer out エラーコード(0:正常終了)
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_PRECONDITION integer in 前処理行列の精度(1:倍精度,2:4倍精度)

5.13 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,iret, CBCG_kstep,CBCG_Eigenflag,power_method_itrmax, CAARNOLDI_sstep,CAARNOLDI_tstep, CAARNOLDI_basis,CAARNOLDI_QRflag )


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

引数名(次元) 入力/出力 説明
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 Compressed Row 形式に格納された行列の非零要素
crsRow_ptr(n+1) integer*4 / integer*8 in Compressed Row 形式のポインタテーブル
crsCol(nnz) integer*4 / integer*8 in Compressed Row 形式に格納された分割行列の非零要素の列番号
ipreflag integer in 前処理指定フラグ(0:無し,1:ポイントヤコビ,2:ブロックヤコビ,3:加法シュワルツ)
ILU_method integer in 不完全LU分解指定フラグ(0:ILU(0),1:D-ILU(対角成分が一致),2:D-ILU(行内の要素和が一致)
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:各プロセス)
iret integer out エラーコード(0:正常終了)
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)

5.14 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,iret, CBCG_kstep,CBCG_Eigenflag,power_method_itrmax, CAARNOLDI_sstep,CAARNOLDI_tstep, CAARNOLDI_basis,CAARNOLDI_QRflag )


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

引数名(次元) 入力/出力 説明
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 Diagonal形式に格納された行列の非零要素
offset(nnd) integer*4 / integer*8 in Diagonal形式に格納された行列の対角からのオフセット
nnd integer*4 / integer*8 in Diagonal形式に格納された行列の対角要素の本数
ipreflag integer in 前処理指定フラグ(0:無し,1:ポイントヤコビ,2:ブロックヤコビ,3:加法シュワルツ)
ILU_method integer in 不完全LU分解指定フラグ(0:ILU(0),1:D-ILU(対角成分が一致),2:D-ILU(行内の要素和が一致)
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:各プロセス)
iret integer out エラーコード(0:正常終了)
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)

5.15 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,iret, CAGMRES_sstep, CAGMRES_tstep, CAGMRES_basis, CAGMRES_QRflag )


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

引数名(次元) 入力/出力 説明
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 Compressed Row 形式に格納された行列の非零要素
crsRow_ptr(n+1) integer*4 / integer*8 in Compressed Row 形式のポインタテーブル
crsCol(nnz) integer*4 / integer*8 in Compressed Row 形式に格納された分割行列の非零要素の列番号
ipreflag integer in 前処理指定フラグ(0:無し,1:ポイントヤコビ,2:ブロックヤコビ,3:加法シュワルツ)
ILU_method integer in 不完全LU分解指定フラグ(0:ILU(0),1:D-ILU(対角成分が一致),2:D-ILU(行内の要素和が一致)
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:各プロセス)
iret integer out エラーコード(0:正常終了)
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)

5.16 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,iret, CAGMRES_sstep, CAGMRES_tstep, CAGMRES_basis, CAGMRES_QRflag )


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

引数名(次元) 入力/出力 説明
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 Diagonal形式に格納された行列の非零要素
offset(nnd) integer*4 / integer*8 in Diagonal形式に格納された行列の対角からのオフセット
nnd integer*4 / integer*8 in Diagonal形式に格納された行列の対角要素の本数
ipreflag integer in 前処理指定フラグ(0:無し,1:ポイントヤコビ,2:ブロックヤコビ,3:加法シュワルツ)
ILU_method integer in 不完全LU分解指定フラグ(0:ILU(0),1:D-ILU(対角成分が一致),2:D-ILU(行内の要素和が一致)
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:各プロセス)
iret integer out エラーコード(0:正常終了)
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)

6 PARCELサブルーチン(間接インターフェース)

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

6.1 PARCEL_KSP_Default_Setting


call PARCEL_KSP_Default_Setting(method)


 PARCEL構造体に初期設定を適用.

引数名(次元) 入力/出力 説明
method type(KSP) out PARCEL 構造体

6.2 set_KSP_CRS


call set_KSP_CRS(icomm,vecx,vecb,n,gn,nnz,istart,crsA,crsRow_ptr,crsCol,itrmax,ipreflag,ILU_method,addL,method)


 Compressed Row形式の行列を設定.

引数名(次元) 入力/出力 説明
icomm integer in MPIコミュニケータ
vecx(n) double precision in 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 Compressed Row 形式に格納された行列の非零要素
crsRow_ptr(n+1) integer*4 / integer*8 in Compressed Row 形式のポインタテーブル
crsCol(nnz) integer*4 / integer*8 in Compressed Row 形式に格納された分割行列の非零要素の列番号
itrmax integer in in:最大反復回数,out:反復回数
ipreflag integer in 前処理指定フラグ(0:無し,1:ポイントヤコビ,2:ブロックヤコビ,3:加法シュワルツ)
ILU_method integer in 不完全LU分解指定フラグ(0:ILU(0),1:D-ILU(対角成分が一致),2:D-ILU(行内の要素和が一致)
addL integer in 加法シュワルツの重なり幅
method type(KSP) in/out PARCEL 構造体

6.3 set_KSP_DIA


call set_KSP_DIA(icomm,vecx,vecb,n,gn,istart,diaA,offset,nnd,itrmax,ipreflag,ILU_method,addL,method)


 Diagonal形式の行列を設定.

引数名(次元) 入力/出力 説明
icomm integer in MPIコミュニケータ
vecx(n) double precision in 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 Diagonal形式に格納された行列の非零要素
offset(nnd) integer*4 / integer*8 in Diagonal形式に格納された行列の対角からのオフセット
nnd integer*4 / integer*8 in Diagonal形式に格納された行列の対角要素の本数
itrmax integer in in:最大反復回数,out:反復回数
ipreflag integer in 前処理指定フラグ(0:無し,1:ポイントヤコビ,2:ブロックヤコビ,3:加法シュワルツ)
ILU_method integer in 不完全LU分解指定フラグ(0:ILU(0),1:D-ILU(対角成分が一致),2:D-ILU(行内の要素和が一致)
addL integer in 加法シュワルツの重なり幅
method type(KSP) in/out PARCEL 構造体

6.4 set_vecx_KSP


 call set_vecx_KSP(n,x,method)


 反復の初期値を設定.

引数名(次元) 入力/出力 説明
n integer*4 / integer*8 in 各プロセスが担当するベクトルの大きさ
x(n) double precision in 反復の初期値
method type(KSP) in/out PARCEL 構造体

6.5 set_vecb_KSP


 call set_vecb_KSP(n,b,method)


 右辺ベクトルbを設定.

引数名(次元) 入力/出力 説明
n integer*4 / integer*8 in 各プロセスが担当するベクトルの大きさ
b(n) double precision in 右辺ベクトル
method type(KSP) in/out PARCEL 構造体

6.6 get_vecx_KSP


 call get_vecx_KSP(n,x,method)


 解ベクトルを取得.

引数名(次元) 入力/出力 説明
n integer*4 / integer*8 in 各プロセスが担当するベクトルの大きさ
x(n) double precision out 解ベクトル
method type(KSP) in PARCEL 構造体

6.7 get_reshistory_KSP


call get_reshistory_KSP(niter,reshistory, method)


 収束履歴を取得する.

引数名(次元) 入力/出力 説明
niter integer out 収束までの反復回数
reshistory(method%itrmax) double precision out 収束履歴
method type(KSP) in PARCEL 構造体

6.8 reset_CRS_Mat


call reset_CRS_Mat(itrmax,ipreflag,ILU_method,addL,crsA,method)


 行列の非零要素位置を変えずにCompressed Row形式の行列値,前処理方法,最大反復回数を変更する.

引数名(次元) 入力/出力 説明
itrmax integer in 最大反復回数(-1:設定済みの値を引き継ぐ)
ipreflag integer in 前処理指定フラグ(-1:設定済みの値を引き継ぐ,0:無し,1:ポイントヤコビ,2:ブロックヤコビ,3:加法シュワルツ)
ILU_method integer in 不完全LU分解指定フラグ(-1:設定済みの値を引き継ぐ,1:対角成分が一致,2:行内の要素和が一致)
addL integer in 加法シュワルツの重なり幅(-1:設定済みの値を引き継ぐ)
crsA(method%nnz) double precision in Compressed Row 形式に格納された行列の非零要素
method type(KSP) in/out PARCEL 構造体

6.9 reset_DIA_Mat


call reset_DIA_Mat(itrmax,ipreflag,ILU_method,addL,diaMat,method)


 行列の非零要素位置を変えずにDiagonal形式の行列値,前処理方法,最大反復回数を変更する.

引数名(次元) 入力/出力 説明
itrmax integer in 最大反復回数(-1:設定済みの値を引き継ぐ)
ipreflag integer in 前処理指定フラグ(-1:設定済みの値を引き継ぐ,0:無し,1:ポイントヤコビ,2:ブロックヤコビ,3:加法シュワルツ)
ILU_method integer in 不完全LU分解指定フラグ(-1:設定済みの値を引き継ぐ,1:対角成分が一致,2:行内の要素和が一致)
addL integer in 加法シュワルツの重なり幅(-1:設定済みの値を引き継ぐ)
diaMat(method%n*method%nnd) double precision in Compressed Row 形式に格納された行列の非零要素
method type(KSP) in/out PARCEL 構造体

6.10 set_KSP_solver


call set_KSP_solver(method,solver)


 クリロフ部分空間法を選択する.

引数名(次元) 入力/出力 説明
method type(KSP) out PARCEL 構造体
solver integer in クリロフ部分空間法フラグ(1:CG法,2:Bi-CGSTAB法,3:GMRES法,4:CAGMRES法,5:CBCG法)

6.11 set_KSP_abstolmax(method,abstolmax)


call set_KSP_abstolmax(method,abstolmax)


 クリロフ部分空間法の絶対残差ノルムによる収束判定条件を設定する.

引数名(次元) 入力/出力 説明
method type(KSP) out PARCEL 構造体
abstolmax double precision in 収束判定条件(絶対残差ノルム)

6.12 set_KSP_rtolmax


call set_KSP_rtolmax(method,rtolmax)


 クリロフ部分空間法の相対残差ノルムによる収束判定条件を設定する.

引数名(次元) 入力/出力 説明
method type(KSP) out PARCEL 構造体
rtolmax double precision in 収束判定条件(相対残差ノルム)

6.13 set_KSP_iovlflag


call set_KSP_iovlflag(method,iovlflag)


 疎行列ベクトル積の通隠蔽手法を選択する.

引数名(次元) 入力/出力 説明
method type(KSP) out PARCEL 構造体
iovlflag integer in 通信隠蔽フラグ(0:無し,1:全プロセス,2:各プロセス)

6.14 set_KSP_GMRES_m(method,GMRES_m)


call set_KSP_GMRES_m(method,GMRES_m)


 GMRES(m)のリスタート数mを選択する.

引数名(次元) 入力/出力 説明
method type(KSP) out PARCEL 構造体
GMRES_m integer in リスタートするまでの反復回数

6.15 set_KSP_GMRES_GSflag


call set_KSP_GMRES_GSflag(method,GMRES_GSflag)


 GMRES(m)の直行化手法を選択する.

引数名(次元) 入力/出力 説明
method type(KSP) out PARCEL 構造体
GMRES_GSflag integer in 直交化フラグ(1:MGS,2:CGS)

6.16 set_KSP_CAGMRES_sstep


call set_KSP_CAGMRES_sstep(method,CAGMRES_sstep)


 CA-GMRES法の省通信ステップ数(s)を選択する.

引数名(次元) 入力/出力 説明
method type(KSP) out PARCEL 構造体
CAGMRES_sstep integer in CA-GMRES法の省通信ステップ数

6.17 set_KSP_CAGMRES_tstep


call set_KSP_CAGMRES_tstep(method,CAGMRES_tstep)


 CA-GMRES法のリスタートステップ数(t)を選択する.

引数名(次元) 入力/出力 説明
method type(KSP) out PARCEL 構造体
CAGMRES_tstep integer in CA-GMRES法のリスタートステップ数

6.18 set_KSP_CAGMRES_basis


call set_KSP_CAGMRES_basis(method,CAGMRES_basis)


 CA-GMRES法の基底生成方法を選択する.

引数名(次元) 入力/出力 説明
method type(KSP) out PARCEL 構造体
CAGMRES_basis integer in CA-GMRES法の基底生成フラグ(0:単基底,1:ニュートン基底)

6.19 set_KSP_CAGMRES_QRflag


call set_KSP_CAGMRES_QRflag(method,CAGMRES_QRflag)


 CA-GMRES法のQR分解方法を選択する.

引数名(次元) 入力/出力 説明
method type(KSP) out PARCEL 構造体
CAGMRES_QRflag integer in CA-GMRES法のQR分解フラグ(1:MGS,2:CGS,3:TSQR,4:CholeskyQR,5:CholeskyQR2)

6.20 set_KSP_CBCG_kstep


call set_KSP_CBCG_kstep(method,CBCG_kstep)


 CBCG法の省通信ステップ数kを選択する.

引数名(次元) 入力/出力 説明
method type(KSP) out PARCEL 構造体
CBCG_kstep integer in CBCG法の省通信ステップ

6.21 set_KSP_CBCG_Eigenflag


call set_KSP_CBCG_Eigenflag(method,CBCG_Eigenflag)


 CBCG法の固有値計算方法を選択する.

引数名(次元) 入力/出力 説明
method type(KSP) out PARCEL 構造体
CBCG_Eigenflag integer in 固有値計算フラグ(1:べき乗法,2:CA-ARNOLDI)

6.22 set_KSP_CAARNOLDI_sstep


call set_KSP_CAARNOLDI_sstep(method,CAARNOLDI_sstep)


 CA-ARNOLDI法の省通信ステップ数(s)を選択する.

引数名(次元) 入力/出力 説明
method type(KSP) out PARCEL 構造体
CAARNOLDI_sstep integer in CA-ARNOLDI法の省通信ステップ数

6.23 set_KSP_CAARNOLDI_tstep


call set_KSP_CAARNOLDI_tstep(method,CAARNOLDI_tstep)


 CA-ARNOLDI法のリスタートステップ(t)数を選択する.

引数名(次元) 入力/出力 説明
method type(KSP) out PARCEL 構造体
CAARNOLDI_tstep integer in CA-ARNOLDI法のリスタートステップ

6.24 set_KSP_CAARNOLDI_basis


call set_KSP_CAARNOLDI_basis(method,CAARNOLDI_basis)


 CA-ARNOLDI法の基底生成方法を選択する.

引数名(次元) 入力/出力 説明
method type(KSP) out PARCEL 構造体
CAARNOLDI_basis integer in CA-ARNOLDI法の基底生成フラグ(0:単基底,1:ニュートン基底)

6.25 set_KSP_CAARNOLDI_QRflag


call set_KSP_CAARNOLDI_QRflag(method,CAARNOLDI_QRflag)


 CA-ARNOLDI法のQR分解を選択する.

引数名(次元) 入力/出力 説明
method type(KSP) out PARCEL 構造体
CAARNOLDI_QRflag integer in CA-ARNOLDI法のQR分解フラグ(1:MGS,2:CGS,3:TSQR,4:CholeskyQR,5:CholeskyQR2)

6.26 set_KSP_power_method_itrmax


call set_KSP_power_method_itrmax(method,power_method_itrmax)


 べき乗法の反復回数を選択する.

引数名(次元) 入力/出力 説明
method type(KSP) out PARCEL 構造体
power_method_itrmax integer in べき乗法の反復回数

6.27 START_KSP_SOLVER


call START_KSP_SOLVER(method)


 連立一次方程式Ax=bの求解.

引数名(次元) 入力/出力 説明
method type(KSP) in/out PARCEL 構造体

7 サブルーチン使用方法

 PARCELのサブルーチン使用方法についてCG法のFORTRANサンプルコードを例にして説明する. Cプログラムにおける利用法の詳細についてはアーカイブのSPARCE/example_Cに格納されているCサンプルコードを参照のこと. Cの場合はFORTRAN用のMPIコミュニケータを用意すればFORTRANと同様の方法でPARCELを利用可能である. CサンプルコードではMPIライブラリの関数MPI_Comm_c2fを用いてCのMPIコミュニケータをFORTRAN用に変換している. 間接呼び出しインターフェースのCサンプルコードでは,FORTRAN側でPARCEL構造体を確保してからC側でそれを制御して利用している.

7.1 Compressed Row形式

 Compressed Row形式のサンプルコードを図8および図9に示す. make_matrix_CRSはPARCELの行列分割に対応したCompressed Row形式の行列を生成する任意のサブルーチンとする. 図8 は直接呼び出しインターフェースを用いた例,図9は間接呼び出しインターフェースを利用した例を示す. 図9ではkrylov,krylov_kernelのモジュールを利用するためSPARSE/srcをコンパイル時にインクルードする必要がある. 非零要素の位置を変更せずに値のみを変更する場合, 図8のコードでは再度parcel_dcgを呼ぶ必要があるため通信相手のリスト作成等の初期化コストが必要となる. 図9のコードではreset_CRS_Matを利用することで非零要素位置の変更がない場合は初期化コストは不要となる. 同様に右辺ベクトル及び反復の初期値を変更する場合も図8ではparcel_dcgを呼ぶため初期化のコストが必要となるが, 図9のコードではset_vecb_KSP,set_vecx_KSPを利用することで初期化コストは不要となる.


  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 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
 
    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 &
         )
    call MPI_Finalize(ierr)
 
    deallocate(vecx)
    deallocate(vecb)
 
    deallocate(reshistory)
 
  end program main


図8 Compressed Row形式 サンプルコード(直接呼び出しインターフェース)

  program main
    use mpi
    use krylov
    use krylov_kernel
    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 iret
    integer ierr
    type(KSP) :: method
 
    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
    vecb=1.0d0
    vecx=1.0d0
 
    call PARCEL_KSP_Default_Setting(method)
    call set_KSP_CRS(&
         MPI_COMM_WORLD, &
         vecx,vecb, &
         n,gn,nnz,istart, &
         crsA,crsRow_ptr,crsCol, &
         itrmax,ipreflag,ILU_method,addL,&
         method)
    call set_KSP_solver(method,1)
    call START_KSP_SOLVER(method)
    call get_vecx_KSP(n,vecx,method)
    call free_KSP(method)
 
    call MPI_Finalize(ierr)
 
    deallocate(vecx)
    deallocate(vecb)
 
    deallocate(reshistory)
 
  end program main


図9 Compressed Row形式 サンプルコード(間接呼び出しインターフェース)

7.2 Diagonal形式

 Diagonal形式のサンプルコードを図10に示す. make_matrix_DIAはPARCELのブロック分割に対応したDiagonal形式の行列を生成する任意のサブルーチンとする. 図10は直接呼び出しインターフェースを用いた例, 図11は間接呼び出しインターフェースを利用した例を示す. 図11の形式ではkrylov,krylov_kernelのモジュールを利用するためSPARSE/srcをコンパイル時にインクルードする必要がある. 非零要素の位置を変更せずに値のみを変更する場合, 図10のコードでは再度parcel_dcg_diaを呼ぶ必要があるため通信相手のリスト作成等の初期化コストが必要となる. 図11のコードではreset_DIA_Matを利用することで非零要素位置の変更がない場合は初期化コストは不要となる. 同様に右辺ベクトル及び反復の初期値を変更する場合も図10ではparcel_dcg_diaを呼ぶため初期化のコストが必要となるが, 図11のコードではset_vecb_KSP,set_vecx_KSPを利用することで初期化コストは不要となる.


  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 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
    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,iret&
         )
 
    call MPI_Finalize(ierr)
 
    deallocate(reshistory)
 
  end program main


図10 Diagonal形式 サンプルコード(直接呼び出しインターフェース)

  program main
    use mpi
    use krylov
    use krylov_kernel
    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 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
    vecb=1.0d0
    vecx=1.0d0
 
    call PARCEL_KSP_Default_Setting(method)
    call set_KSP_DIA(&
         MPI_COMM_WORLD, &
         vecx,vecb, &
         n,gn,istart, &
         diaA,offset,nnd,&
         itrmax,ipreflag,ILU_method,addL,&
         method)
    call set_KSP_solver(method,1)
    call START_KSP_SOLVER(method)
    call get_vecx_KSP(n,vecx,method)
    call free_KSP(method)
 
    call MPI_Finalize(ierr)
 
    deallocate(vecx)
    deallocate(vecb)
 
    deallocate(reshistory)
 
  end program main


図11 Diagonal形式 サンプルコード(間接呼び出しインターフェース)