筑波大学計算科学研究センター
CCS HPCサマーセミナー「並列数値アルゴリズム I」
多田野 寛人tadano@cs.tsukuba.ac.jp

筑波大学計算科学研究センター

講義内容
・連立一次方程式の求解法- 様々な行列に対する解法- 疎行列の扱い方- 基本的な線形演算の MPI での並列化

・複数本の右辺ベクトルをもつ連立一次方程式の解法- 解法の特色- 効率的な計算手法と OpenMP での並列化

・レポート課題について-1-

連立一次方程式の求解法-2-

自然現象・工学現象の解析
自然現象・工学現象モデル化
偏微分方程式の初期値・境界値問題

偏微分方程式の近似解
解析

離散化

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

連立一次方程式は様々な分野における数値シミュレーションで現れ，計算時間の大部分が求解に費やされている-3-

連立一次方程式
連立一次方程式：Ax = b2666 a1,1666 a
6 2,1A = 66666 .
666 ..
4an,1a1,2a2,2..
.

...
...
..
.

an,2...

3232
3a1,n 77xb6
777
7611666
777666777
666 x2 777
666 b2 7777
a2,n 7767776.. 7777 , x = 6666 .. 7777 , b = 6666 .. 7777
666 . 777666 . 777
. 77754554
an,nxnbn

連立一次方程式は様々な分野における数値シミュレーションで現れ，計算時間の大部分が求解に費やされている-4-

直接解法と反復解法
直接解法• Gaussの消去法，LU分解法など• 有限回の演算で必ず解くことができる• 係数行列 A を変形するため，非零要素数が増大

反復解法• 定常反復法，クリロフ部分空間反復法• 必要な演算は係数行列 A とベクトルの積，内積など
係数行列の疎性をそのまま使える• 問題によっては多くの反復回数を要することがある-5-代表的な直接解法（Gaussの消去法）
Gaussの消去法では，連立一次方程式 Ax = b を変形し，上三角行列 U をもつ連立一次方程式にした後で解ベクトル x を計算する．
2666 a1,1666 a666 2,1666 ..
666 .
4an,1

連立一次方程式： Ax = ba1,2a2,2..
.

...
...
..
.

an,2...

32
3 2
3a1,n 77 66 x1 77 66 b1 77
76
77 6a2,n 7777 6666 x2 7777 6666 b2 7777
77 676.. 7777 6666 .. 7777 = 6666 .. 7777
. 7775 6664 . 7775 6664 . 7775
xnan,nbn

変形

2666 u1,1666666
666
666
40

変形後の方程式： U x = b0u1,2u2,2...
...
..
.

3233 2u1,n 77 66 x1 77 66 b01 77
76
7 67u2,n 7777 6666 x2 7777 6666 b02 7777
76
7 67.. 7777 6666 .. 7777 = 6666 .. 7777
. 7775 6664 . 7775 6664 . 7775b0nxnun,n変形後は下記の後退代入計算により，解ベクトル x の成分を計算する．
⌘1 ⇣ 0b − u j, j+1 x j+1 − u j, j+2 x j+2 − · · · − u j,n xn , j = n, n − 1, . . . , 1.
xj =u j, j j

上三角行列をもつ方程式への変形過程の計算量は O(n3 ) ，後退代入の計算量は O(n2 ) ．
-6-代表的な直接解法（LU分解法）
LU分解法は係数行列 A を下三角行列 L と上三角行列 U の積に分解することで，Ax = b を解く方法である．係数行列だけが変形され，右辺ベクトルは変わらない．
2666 1666 l666 2,1666 ..
666 .
4ln1

変形後の方程式： LU x = b1..
.
ln,2解ベクトルの計算手順

..

.
...

0132
777 666 u1,1777 666
777 666
777 666
777 666
540u1,2u2,2...
...
..
.

32
3 23u1,n 77 66 x1 77 66 b1 77
76
7 67u2,n 7777 6666 x2 7777 6666 b2 7777
767 67.. 7777 6666 .. 7777 = 6666 .. 7777
. 7775 6664 . 7775 6664 . 7775
xnbnunn

1. y = U x とおいて Ly = b を前進代入で解き，y を求める．
2. U x = y を後退代入で解く．

行列 A のLU分解と解ベクトルの計算量はそれぞれ，O(n ), O(n ) である．
右辺ベクトルが複数ある場合は，LU分解法を用いた方が効率良く求解できる．
32-7-疎行列の直接解法ソフトウェア– MUMPS（MPI）• http://mumps.enseeiht.fr/– SuperLU（SuperLU_DIST で MPI並列も可）• http://crd-legacy.lbl.gov/~xiaoye/SuperLU/– Pardiso（MPI, スレッド並列も可）• Intel Math Kernel Library (MKL) にも入っている• http://www.pardiso-project.org/– WSMP（MPI, スレッド並列も可）• http://researcher.watson.ibm.com/researcher/view_project.php?id=1426-8-クリロフ部分空間とは？
クリロフ部分空間

行列 A と非ゼロベクトル v から生成されるベクトル列で張られる部分空間：
Kk (A, v) ≡ span v, Av, . . . , A

k−1


v

をクリロフ部分空間という．

クリロフ部分空間反復法

アレクセイ・クリロフ（Wikipediaより）• クリロフ部分空間を用いて連立一次方程式 Ax = b の近似解を生成する方法をクリロフ部分空間反復法という．
• 演算の主要部は行列 A とベクトルの積で，計算量は αn2 .
（ α は係数行列の非零要素数の割合）-9-クリロフ部分空間反復法• 任意の初期解 x0 を与え，反復によって近似解を更新し，真の解 x を探索• r0 は初期残差と呼ばれ，r0 = b – Ax0 で計算される• 一般に真の解 x は分からないため，残差 rk = b – Axk のノルムを用いて真の解への収束をチェックする

x0 + Kk (A; r0 )x0 + Kk−1 (A; r0 )x0 + K3 (A; r0 )x0 + K2 (A; r0 )x0 + K1 (A; r0 )•x1x0•

•

x2...

• x3

クリロフ部分空間反復法の概念図

xk

•
•

x

•

xk−1

- 10 -

エルミート行列に対する解法
1. 係数行列がエルミート行列（A = AH）の場合
・共役勾配法（Conjugate Gradient method: CG法）
・共役残差法（Conjugate Residual method: CR法）

・最小残差法（Minimal Residual Method: MINRES法）
係数行列のエルミート性を使うことで短い漸化式（計算量が少ない）のアルゴリズムが導出できる

補足：エルミート行列
A = AH = ĀT(ai j = ā ji )- 11 -

共役勾配法（CG法）のアルゴリズム
x0 is an initial guess,
Compute r0 = b − Ax0 ,Set p0 = r0 ,For k = 0, 1, . . . , until rk 2 ≤ εTOL b2 do :qk = Apk ,行列ベクトル積(rk , rk ),αk =内積( pk , qk )xk+1 = xk + αk pk ,ベクトルの定数倍と加算
rk+1 = rk − αk qk ,(rk+1 , rk+1 ),βk =内積(rk , rk )pk+1 = rk+1 + βk pk ,ベクトルの定数倍と加算
End For- 12 -

相対残差, rk 2 /b2

CG法の相対残差履歴のイメージこのイメージでは
rk 2 /b2 ≤ 10−12を満たしたら停止反復回数, k・反復の過程で，相対残差（ rk 2 /b2 ）をチェックする．
・自分で定めた条件を満たしたら反復をやめ，xk を解として採用．

- 13 -

非エルミート行列に対する解法
2. 係数行列が 非エルミート行列（ A  AH ）の場合
短い漸化式を用いる解法
・BiCG法・CGS法・BiCGSTAB法
・GPBiCG法 ・BiCGSTAB(  )法 ・IDR(s)法など

計算量は少ないが，残差ノルムは単調減少しない

長い漸化式を用いる解法
・GCR法

・GMRES法

残差ノルムは単調減少するが，計算量が多い- 14 -

BiCG法のアルゴリズム

x0 is an initial guess,
Compute r0 = b − Ax0 ,Choose r∗0 such that (r∗0 , r0 )  0,Set p0 = r0 and p∗0 = r∗0 ,For k = 0, 1, . . . , until rk 2 ≤ εTOL b2 do:q∗k = AH p∗k ,qk = Apk ,(r∗k , rk ),αk = ∗( pk , qk )xk+1 = xk + αk pk ,rk+1 = rk − αk qk , r∗k+1 = r∗k − ᾱk q∗k ,(r∗k+1 , rk+1 ),βk =(r∗k , rk )pk+1 = rk+1 + βk pk , p∗k+1 = r∗k+1 + β̄k p∗k ,End For

行列ベクトル積
内積
ベクトルの定数倍と加算- 15 -BiCGSTAB法のアルゴリズム
x0 is an initial guess,
Compute r0 = b − Ax0 ,Choose r∗0 such that (r∗0 , r0 )  0,Set p0 = r0 ,For k = 0, 1, . . . , until rk 2 ≤ b2 do:
qk = Apk ,(r∗0 , rk ),αk = ∗(r0 , qk )t k = rk − αk qk ,sk = At k ,(sk , t k )ζk =,(sk , sk )xk+1 = xk + αk pk + ζk t k ,rk+1 = t k − ζk sk ,αk (r∗0 , rk+1 ),·βk =
ζk (r∗0 , rk )pk+1 = rk+1 + βk ( pk − ζk qk ),End For

行列ベクトル積
内積ベクトルの定数倍と加算- 16 -

GCR法のアルゴリズム

x0 is an initial guess,
Compute r0 = b − Ax0 ,Set p0 = r0 and q0 = s0 = Ar0 ,For k = 0, 1, . . . , until rk 2 ≤ εTOL b2 do :(qk , rk ),αk =(qk , qk )xk+1 = xk + αk pk ,rk+1 = rk − αk qk ,sk+1 = Ark+1 ,(qi , sk+1 )βk,i = −, (i = 0, 1, . . . , k)(qi , qi )kpk+1 = rk+1 + ∑ βk,i pi ,i=0kqk+1 = sk+1 + ∑ βk,i qi ,End For

i=0・行列ベクトル積は1反復あたり1回・長い漸化式を使うため多くのメモリが必要
・リスタートにより演算量とメモリ量を削減- 17 -

反復法の収束特性
相対残差, rk 2 /b2

行列 A：poisson3Da（n = 13,514，nnz(A) = 352,762, SuiteSparse Matrix Collection）
右辺 b ：b = [1, 1, . . . , 1]T

反復回数, k反復法の相対残差履歴■：BiCG法，■：CGS法，■：BiCGSTAB法，■：GMRES法- 18 -

複素対称行列に対する解法
3. 係数行列が複素対称行列（ A = AT  AH ）の場合
共役直交共役勾配法（Conjugate Orthogonal Conjugate Gradient: COCG法）

係数行列が複素対称行列の場合は，1反復あたり1回の行列ベクトル積で，かつ短い漸化式で計算ができる

補足：複素対称行列
A = AT  AH(ai j = a ji  ā ji )- 19 -

COCG法のアルゴリズム
x0 is an initial guess,Compute r0 = b − Ax0 ,Set p0 = r0 ,For k = 0, 1, . . . , until rk 2 ≤ εTOL b2 do :qk = Apk ,行列ベクトル積(r̄k , rk )
,
αk =
内積
( p̄k , qk )xk+1 = xk + αk pk ,ベクトルの定数倍と加算
rk+1 = rk − αk qk ,(r̄k+1 , rk+1 )
,
βk =
内積
(r̄k , rk )pk+1 = rk+1 + βk pk ,ベクトルの定数倍と加算
End For- 20 -

疎行列が現れる例
2次元 Poisson 問題 22∂uu∂+= f , in Ω22∂x∂y u = ū,on ∂Ωf, ū は既知の関数．
W をx, y 方向に (M+1) 等分し5点中心差分で離散化

y1

O1x

M×M次行列をもつ連立一次方程式に帰着
行列の要素数：M4 個
非零要素の数：5M2 - 4M 個

- 21 -

疎行列の格納形式（CRS形式）
Compressed Row Storage (CRS) 形式【補足】Compressed Sparse Row (CSR) と呼ぶこともあります．
val: 非零要素を格納する配列 a11 0 a13 0 a15  0 a22 0 a24 a25 A =  a31 a32 a33 0 0  0 0 a a
43 44 0 
0 a52 0 a54 a55

col_ind: 非零要素の列番号を格納
row_ptr: 各行の先頭の非零要素が格納されている場所を指す配列

val: a11 a13 a15 a22 a24 a25 a31 a32 a33 a43 a44 a52 a54 a55
col_ind:13524

row_ptr:14710125

1233424515

最後は非零要素数+1 の値- 22 -疎行列の格納形式（CCS形式）
Compressed Column Storage (CCS) 形式【補足】Compressed Sparse Column (CSC) と呼ぶこともあります．
val: 非零要素を格納する配列

 a11 0 a13 0 a15  0 a
22 0 a24 a25 
A =  a31 a32 a33 0 0  0 0 a a
43 44 0 
0 a52 0 a54 a55

row_ind: 非零要素の行番号を格納
col_ptr: 各列の先頭の非零要素が格納されている場所を指す配列

val: a11 a31 a22 a32 a52 a13 a33 a43 a24 a44 a54 a15 a25 a55
row_ind:132351

col_ptr:

1369121534245125

最後は非零要素数+1 の値- 23 -CRS形式の行列ベクトル積
CRS形式では，行列 A の非零要素は行方向に連続にアクセスされる．即ち，y の i 番目の成分 yi は，行列 A の i 行目とベクトル x の内積で求められる．

行列 A とベクトル x の積 y = Ax  y1   a11 a12 . . . a1n   x1  y   a 2   21 a22 . . . a2n   x2 
 ..  =  ..
..
..   .. ..
 .   .
  . .
.
.
 
yn
an1 an2 . . . ann
xn

Fortran のコード例
do i=1,ny(i) = 0.0D0
do j=row_ptr(i), row_ptr(i+1)-1y(i) = y(i)+val(j)*x(col_ind(j))
end doend do- 24 -

CCS形式の行列ベクトル積
CCS形式では，行列 A の非零要素は列方向に連続にアクセスされる．即ち，行列ベクトル積 y = Ax は，行列 A の列ベクトルの線型結合として求められる．
行列 A とベクトル x の積 y = Ax x1  x n 2 y = [a1 , a2 , . . . , an ]  .  = ∑ ai xi
 ..  i=1
xn

Fortran のコード例
do i=1,ny(i) = 0.0D0
end dodo j=1,ndo i=col_ptr(j),col_ptr(j+1)-1k = row_ind(i)y(k) = y(k) + val(i) * x(j)end doend do- 25 -

CRS形式の転置行列ベクトル積
CRS形式の転置行列ベクトル積 z = AT x2666 a1,1666 0666z = 6666 a1,3666 064a1,5|0a2,20a2,4a2,5a3,1a3,2a3,300{zAT

00a4,3a4,4032
30 77 66 x1 77
76
7a5,2 7777 6666 x2 7777
76
70 7777 6666 x3 7777
76
7a5,4 77775 66664 x4 77775
a5,5x5} | {z }xFortran のコード例
do i=1,nz(i) = 0.0D0
end dodo j=1,ndo i=row_ptr(j),row_ptr(j+1)-1
k = col_ind(j)z(k) = z(k) + val(i) * x(i)end doend do• CRS形式で格納される行列の非零要素は，行方向に連続にアクセスされる．
• この行列の転置をとると，非零要素は列方向に連絡にアクセスされることになる．
• 即ち，CRS形式で格納された行列の転置行列ベクトル積を行うには，CCS形式の行列ベクトル積コードを実行すればよい，ということになる．
- 26 -

CCS形式の転置行列ベクトル積
CCS形式の転置行列ベクトル積 z = AT x2666 a1,1666 0666z = 6666 a1,3666 0
64a1,5|0a2,20a2,4a2,5a3,1a3,2a3,300{zAT

00a4,3a4,4032
30 77 66 x1 77
76
7a5,2 7777 6666 x2 7777
76
70 7777 6666 x3 7777
76
7a5,4 77775 66664 x4 77775
a5,5x5} | {z }Fortran のコード例
do i=1,nz(i) = 0.0D0
do j=row_ptr(i), row_ptr(i+1)-1z(i) = z(i)+val(j)*x(col_ind(j))
end doend dox• CRS形式では行列の非零要素は列方向に連続にアクセスされる，転置をとると非零要素は行方向に連絡にアクセスされることになる．
• 先程と同様に，CCS形式で格納された行列の転置行列ベクトル積を行うには，CRS形式の行列ベクトル積コードを実行すればよい．
- 27 -

行列ベクトル積の並列化
・CRS形式の場合の y = Ax の計算
Proc. 0
Proc. 1
Proc. 2=*

Proc. 3

Axy

各プロセスで分散してデータを保持する全てのプロセスで保持する

MPI_Gather でProc. 0 に集める
- 28 -

行列ベクトル積の並列化

Proc. 3

Proc. 2

Proc. 1

Proc. 0・CCS形式の場合の y = Ax の計算=*+++Axy

各プロセスで分散してデータを保持する

分散して保持する

MPI_Reduce でProc. 0に結果を足し合わせて送る- 29 -

Proc. 1=j=1

x jy jProc. 2

Proc. 3=

Proc. 0=x(x, y) =n=

内積の並列化ytmp_sum tmp_sum

tmp_sum tmp_sum

sum
MPI_Reduce で Proc.0 に集める

- 30 -

内積計算の MPI コード例

n
program main(x, y) =x jy jinclude 'mpif.h'j=1...
call mpi_init(ierr)call mpi_comm_size(mpi_comm_world, nprocs, ierr)call mpi_comm_rank(mpi_comm_world, myrank, ierr)...
tmp_sum = 0.0D0
do i=istart(myrank+1), iend(myrank+1)
tmp_sum = tmp_sum + x(i) * y(i)end docall mpi_reduce(tmp_sum, sum, 1, mpi_double_precision,mpi_sum, 0, mpi_comm_world, ierr)...
call mpi_finalize(ierr)- 31 -

ベクトルの定数倍と加算の並列化
y = y + ax （x, y：ベクトル，a：スカラー）a を MPI_Bcast で全プロセスに送る
Proc. 0a

Proc. 1

Proc. 2

Proc. 3

aaaa++++xy- 32 -

複数本の右辺ベクトルをもつ連立一次方程式の解法- 33 -

複数右辺ベクトルをもつ方程式
右辺が s 本の連立一次方程式
AX = B

ここで，A : n次行列，
(1) (2)(s)(1) (2)(s)
X = x ,x ,...,x ,B = b ,b ,...,b
直接法による解法
・係数行列の完全分解（A = LU など）が必要
・完全分解できれば，s 回の前進・後退代入で OK・分解には多くの計算量，メモリ量が必要
反復法で s 本の方程式を効率よく解けないか？
- 34 -

ブロッククリロフ部分空間反復法
ブロッククリロフ部分空間反復法の種類

・Block BiCG法
・Block GMRES法・Block QMR法・Block BiCGSTAB法
・Block BiCGGR法

O’Leary (1980)
Vital (1990)
Freund (1997)
El Guennouni (2003)
Tadano (2009)

複数本の右辺ベクトルをまとめて扱うことで効率よく解を求めることができる- 35 -

ブロッククリロフ部分空間反復法「効率よく」とはどういうことか？

相対残差, Rk F /BF

1本ずつ解いた場合より少ない反復回数で残差が収束する（ことがある）

反復回数, kBlock BiCGSTAB法の相対残差履歴．
ここで，■：s = 1，■：s = 2，■：s = 4．

- 36 -

Block CG法
X0 ∈ Cn×s is an initial guess,
Compute R0 = B − AX0 ,Set P0 = R0 ,For k = 0, 1, . . . , until Rk F ≤ εTOL BF do:
Qk = APk ,H
Solve (PHk Qk )αk = Rk Rk for αk ,Xk+1 = Xk + Pk αk ,Rk+1 = Rk − Qk αk ,H
Solve (RHk Rk )βk = Rk+1 Rk+1 for βk ,CG法と異なるところ1. 行列・ベクトル積の本数が1本から s 本に増加．
2. αk , βk が s 次行列になった．
3. ベクトルの定数倍の計算が行列・行列積になった．

Pk+1 = Rk+1 + Pk βk ,End For- 37 -

行列・ベクトル積の効率化

・行列は CRS 形式で格納されているとする．
・Y = AX を計算．Y, X は n 行 s 列の配列．
do k=1,s! s についてのループが一番外側
do i=1,ndo j=row_ptr(i), row_ptr(i+1)-1Y(i,k) = Y(i,k) + val(j) * X(col_ind(j),k)end doend doend do[ 問題点 ]・X についてメモリの連続アクセスができていない（Fortran は行方向に連続にデータが並んでいる）
・行列データを s 回読まなければならない．

- 38 -行列・ベクトル積の効率化[ 解決策 ]・X, Y を転置した形で保持させる
do i=1,ndo j=row_ptr(i), row_ptr(i+1)-1
do k=1,s ! s についてのループが一番内側
Y(k,i)= Y (k,i) + val(j) * X(k,col_ind(j))
end doend doend do・X について（少なくとも s 回は）連続アクセスができる．
・行列データは1回しか読み込まない．
・Y についても連続アクセスが保たれている．

- 39 -

n×s 行列・s×s 行列積の計算

・行列・ベクトル積の効率化のためにベクトルを転置した．
・n×s 行列と s×s 行列の積の計算も工夫が必要．
Xk+1 = Xk + Pk αk

転置T
Xk+1= XkT + αTk PTk

do j=1,ndo i=1,sdo k=1,sX(k,j) = X(k,j) + Alpha(k,i) * P(i,j)
end doend doend doAlpha は転置済みとする．

ベクトルと αk を全て転置して保持することで連続アクセスが可能に．

- 40 -

s×n 行列・n×s 行列積の計算
・ak, bk を求めるために必要．

・ Ck = PHk Qk を計算することを考える．
do j=1,ndo i=1,s
do k=1,sC(k,i) = C(k,i) + dconjg(P(k,j)) * Q(i,j)
end doend doend do- Ck の計算も，連続メモリアクセスを保持できる．
- dconjg：複素数の共役を求める関数．
- 41 -

OpenMP による並列化
・共有メモリ向けの並列化インターフェース．
・既存のプログラムに数行加えるだけで並列化ができる．
!$OMP PARALLEL【プログラム】!$OMP END PARALLELと書くと，スレッドが立ち上がり，スレッド毎に別々の処理ができるようになる．
（以下のコードは，!$OMP PARALLEL と!$OMP END PARALLEL で囲まれているとします．）- 42 -

OpenMP による並列化
1. 行列・ベクトル積の並列化!$OMP DO PRIVATE(j,k)
do i=1,ndo j=row_ptr(i), row_ptr(i+1)-1
do k=1,sY(k,i)=Y(k,i)+val(j)*X(k,col_ind(j))
end doend doend do!$OMP END DO

最初の do ループの前に !$OMP DO ... を加えるだけ．
最後の !$OMP END DO は省略可能．

- 43 -

OpenMP による並列化
2. n×L行列・L×L行列の積の並列化!$OMP DO PRIVATE(i,k)
do j=1,ndo i=1,sdo k=1,sX(k,j) = X(k,j) + Alpha(k,i) * P(i,j)
end doend doend do!$OMP END DO

一番外の do ループの前に !$OMP DO ... を加えるだけ．
最後の !$OMP END DO は省略可能．

- 44 -

OpenMP による並列化

3. L×n行列と n×L行列の積の並列化
以下を実行して Ck = PHk Qk を計算する．

!$OMP SINGLE
do j=1,L
do i=1,L
C(i,j) = 0.0D0
end doend do!$OMP END SINGLE!$OMP DO PRIVATE(i,k) REDUCTION(+:C)
do j=1,ndo i=1,s
do k=1,sC(k,i) = C(k,i) + dconjg(P(k,j)) * Q(i,j)
end doend doend do!$OMP END DO- 45 -

実行時間 / s [秒]行列ベクトル積の性能改良版
単純な実装

右辺ベクトル数, s右辺ベクトル数の変化に対する計算時間変化（OpenMP 12並列）．
・行列：格子量子色力学計算（QCD）で現れる行列
・行列サイズ：1,572,864，非ゼロ要素数：80,216,064．
・実験環境：CPU : Intel Xeon E5-2620v3 2.4GHz × 2，・コンパイラ：gfortran ver. 5.4，オプション：-O3 –fopenmp- 46 -

まとめ
・連立一次方程式の解法であるクリロフ部分空間反復法を取り上げた．

・疎行列に対する行列・ベクトル積の実装
方法とその並列化について述べた．

・ブロッククリロフ部分空間反復法とコードの最適化，及び OpenMP での並列化について述べた．
- 47 -

レポート課題

次の行列 A と右辺ベクトル b をもつ連立一次方程式 Ax = b を，BiCG法とBiCGSTAB法で解くプログラムを作りなさい．但し，行列 A はCRS形式で格納すること．
行列サイズ n は任意とする．
 2 γA = 



012..
.

1..
.
γ..

.
2γ0
 
 1 
 1  
 . 
 , b = A  .. 
 1 
 1 12−12• 行列 A のパラメータ は 0 < < 1 とする．反復は条件：krk k2 /kbk2  10 を満たしたら停止すること．初期解 x0 はゼロベクトルとする．
• 複数のパラメータ について実験し，反復過程における相対残差 krk k2 /kbk2 をグラフにプロットするとともに考察を述べなさい．
• プログラムリストも提出すること．プログラミング言語は何を用いてもよい．

- 48 -