JPH0666069B2 - ベクトル計算機の記憶装置への大規模不規則疎行列格納法 - Google Patents

Description

【発明の詳細な説明】 （産業上の利用分野） 本発明は数値計算、特に有限要素法プログラムにおける
大規模不規則疎行列の格納法に関し、特に、CG法・チェ
ビシェフ反復法等行列とベクトルの積算が主たるコスト
を占める解法を用いるときの、行列各要素の、ベクトル
計算機の記憶装置上への大規模不規則疎行列格納法に関
するものである。

（従来の技術） 従来は、第３図（ａ）〜（ｆ）のように、行列の非ゼロ
要素を行毎にまとめて一次元化する格納法が使われてき
た。これは、行列Ａが第３図（ａ）の非ゼロ／ゼロ要素
のパターンを持つとするとき、行列の非ゼロ要素（ｘ）
のみを、第３図（ｂ）の→に沿って１行目から行毎に、
配列Ａに格納する方法である。行列Ａが対称行列のとき
は第３図（ｃ）のように上三角部のみを格納する。第３
図（ｄ）は、このようにして格納したときの、配列Ａの
先頭の数要素を表している。配列の各要素と、元の行列
との対応を知るために、別に、整数配列JA、IAを使う。
行列のｉ行目の非ゼロ要素の先頭の番地をIA（ｉ）に格
納する。配列ＡのIA（ｉ）番地からIA（ｉ＋１）−１番
地ｉ行目の非ゼロ要素が対応することになる。配列JAの
各要素は、同一番地のＡの要素の元の行列での列番号を
示す。第３図（ｅ）及び（ｆ）は、それぞれ、IA及びJA
の先頭の部分を表している。この格納法を用いるとき
の、行列×ベクトルのFortranコードが第４図である。
最深DO LOOPが各ｉについて、 を実行する。

さらにベクトル計算機の特性を考慮した格納法として、
第５図のような方法が用いられている。第５図（ａ）
は、格納しようとする行列Ａの非ゼロ／ゼロパターン、
第５図（ｂ）及び（ｃ）は、その上三角部及び下三角部
である。この格納法では、行列の次元数をＮ、下三角部
の一行当たりの非ゼロ要素数の最大値をMLとして、Ｎ×
MLの２次元配列ALを用意し、第５図（ｂ）の下三角部の
ゼロ要素を左方向に圧縮して、２次元配列ALに格納す
る。各行に於て、要素数がMLに満たない場合は、行の後
ろにゼロ要素を詰める（第５図（ｄ）参照）。また、元
の行列での列番号の確保の為に、同じ大きさの２次元配
列JLを用意して、各要素の列番号にこれを格納する。ゼ
ロ要素を詰めたところは、適当な列番号を詰める。

同様に、NU×MUの２次元配列AU、JUを用意し、上三角部
の最後のゼロ要素だけからなる行を除いたものを左詰め
に格納し、対応する列番号をJUに格納する（第５図
（ｅ）参照）。ここで、MUは上三角部の一行当たりの非
ゼロ要素の最大値、NUはＮから最後の、ゼロ行の数を除
いたものである。この格納法のための行列×ベクトルの
Fortranコードが第６図で、（DO BLOCK 10）で行列の対
角要素の寄与を計算し、（DO BLOCK 20）で上三角部の
寄与を、（DO BLOCK 30）で下三角部の寄与を計算して
いる。

（発明が解決しようとする課題） 従来方法である第３図の方法では、最深ループ長が十分
でない。また第５図の方法では、ループ長を延ばす為に
余分なゼロ要素を数多く詰めており、特に一行当たりの
非ゼロ要素数の平均がその最大値よりかなり小さい場合
には、計算量が過大になりベクトル化の効果が期待でき
ない。また、同一列中に同じ列添字が並ぶことが多いた
め、メモリアクセス競合が発生しやすい。以上の理由
で、従来の方法ではベクトル計算機の特性が十分には引
き出されておらず、計算コストの増加を招く。

（課題を解決するための手段） 本発明のベクトル計算機の記憶装置への大規模不規則疎
行列格納法では、ベクトル長を十分長くとるため、行列
の非ゼロ要素を左方に圧縮して詰め、行列・ベクトル積
の計算時に列方向にアクセスする。この際、余分のゼロ
要素を詰める必要が生じないよう、予めベクトル要素の
番号付け方を、行毎の非ゼロ要素が単調減少になるよう
に付け替えている。また、バンクコンフリクトが起きな
いよう、行列を上三角、対角、下三角の三部に分け、対
角部分だけは別に格納し、下三角と上三角について、上
三角、下三角の順に行毎に左詰めし、縦方向に一次元化
している。

（実施例） 以下本発明の実施例において図面を参照しながら説明す
る。第１図（ａ）〜（ｆ）は、本発明の格納法を具体的
に表した図で、（ａ）は、格納しようとする行列Ａのゼ
ロ／非ゼロパターンを示している。ここで、各行の非ゼ
ロ要素の数は降順になっているものとする。第１図
（ｂ）及び（ｃ）は、その上三角部及び下三角部であ
る。以下、格納法の説明のために、行列と同一サイズの
２次元配列を考え、その配列に、上三角部を左に圧縮し
て詰め（第１図（ｄ））、次に、下三角部を左に圧縮す
る（第１図（ｅ））。次に、第１図（ｆ）の↓に沿って
１列目から、順に配列Ａに、非ゼロ要素を、１次元的に
格納していく。各要素の列番号は別の整数配列JAに順に
格納して行く。また、第１図（ｅ）の各列の先頭の要素
の番地を配列JHに格納する。対角要素は別の配列に格納
する。

第２図は、この格納法での行列×ベクトルのFortranコ
ードで、（DO BLOCK 10）で、対角要素の寄与を計算
し、（DO BLOCK 20）は各ｊについて、第１図（ｅ）で
の各ｊ列目の要素の寄与を計算している。

第８図は、第３図（ｃ）の従来の格納法から、本格納法
に変換するアルゴリズムのフローである。Ｎは、ベクト
ルの長さを、Ndは、行列の非ゼロ要素数を表す。IAは第
３図（ｃ）の各行の先頭の番地を保持している。また、
過程1.1、２、３は、Ｎに、過程1.2、４は、Ndに比例す
る演算量だけを含んでいる。とくに、過程４について
は、 よりこのことがわかる。従って、全体で、Ndに比例する
演算量以下で変換が可能である。

次に、本実施例の全体図に移る。第９図は、単位円で
の、ラプラシアンの最小固有値問題を有限要素法で離散
化し、行列の固有値問題に還元して求める際に本格納法
を用いた例のフローである。計算の手順は、 （１）領域の分割を行う。本実施例では、第10図のよう
に分割している。

（２）初期節点番号付けを行う。

（３）それを、行列を作ったときの行毎の非ゼロ要素数
が降順になるように節点番号を付け直す。先ず、各行の
非ゼロ要素数を計算し、それが行列の次元数Ｎに比べて
小さいことを利用し、第７図の線形ソート法を用いて付
け替える。このアルゴリズムの採用によって、節点番号
の付け替えは、ベクトルの長さＮの数倍程度の演算回路
で可能になる。

（４）行方向格納時の、各行の先頭の番地、列番号をも
とめる。本実施例では、生成される行列が対称のため、
第３図（ｃ）の対称行列用行方向格納法を用いている。

（５）従来の行毎の格納法第３図（ｃ）による係数行列
生成を行う。

（６）それを第８図の変換の手段を用いて、第１図の本
発明の格納法に変換する。

（７）このあと、CG法（Conjugate Gradient Method）
を用いてラプラシアンの最小固有値を求めて、計算を終
了する。CG法では、一反復あたり、４回の行列・ベクト
ルの積、７回のベクトルの内積を多数回繰り返す。本実
施例では、比較のため、行列×ベクトルを実行する際
に、第１図の本発明によるものと、前述の、第３図及び
第５図の格納法によるもので行っている。

本実施例７のCG法に於て現れる行列×ベクトルの演算時
間の、１非ゼロ要素当たりの時間を各行列の格納法別に
プロットしたのが第11図である。縦軸に１非ゼロ要素当
たりの計算時間、横軸に有限要素法の基底関数の次数を
とっている。図中に示されているT1Pが本格納法の結
果、Ｒが行毎の従来法、Ｈが第５図の方法である。基底
関数の次数が低い場合、即ち１行中に非ゼロ要素数がよ
り少ない場合、行方向格納法は非常にコストが高い。一
方、次数が高い場合は、即ち１行中の非ゼロ要素数にば
らつきがある場合は、第５図の方法は、コストが高くな
る。本格納法によるものは、何れの次数でも、従来の方
法より速く、且つ、安定した値を示しており、平均的
に、行毎のものの10倍、第５図の方法の８倍程度速くな
っている。

第12図は、第９図のフローで与えられる、本実施例全体
の実行時間を各行列格納法別に計測したもので、左から
順に、節点再番号付け及び格納法変換のコストを含む第
９図過程１〜６の行列作成のコスト、第９図過程７のCG
方のコスト、それに、全体のコストを各限要素法の次数
１〜４毎に表にしたものである。CG法のコスト低下は、
１回の行列・ベクトル積のコストの低減×（４×CGB法
の反復回数＋２）によって、第11図と関連する。本発明
の格納法を用いるために要する、行列作成のコストの増
加は、本発明の格納法を用いることによるCG法のコスト
の減少に比べて、僅かであり、全体のコストは、従来の
ものの、２〜７割程度に減少している。

（発明の効果） 以上説明したように、本発明のベクトル計算機上でのた
とえば乗算向き大規模疎行列格納法では、行毎の要素を
左詰めにし列毎にアクセスするため、１行当たりの非ゼ
ロ要素数が小さいときでも、最深ループのループ長が長
い。このとき、予め一行中の非ゼロ要素が降順になるよ
うに節点番号をつけ替えてあるために、ゼロ要素を詰め
る必要がなく、１行当たりの非ゼロ要素数の最大値が、
その平均値に比べて大きいときも、計算量の増加をもた
らさない。また、行列を上三角、下三角の順で左詰めす
ることにより、各列同一要素を参照することが希になり
メモリアクセス競合を防ぐことができる。

【図面の簡単な説明】

第１図は本発明のベクトル計算機の記憶装置への不規則
疎行列格納法を示す図で、（ａ）は行列Ａのゼロ・非ゼ
ロパターンを示す図、（ｂ）は同上三角部のみを示す
図、（ｃ）は同下三角部のみを示す図、（ｄ）は上三角
部を左に圧縮した図、（ｅ）は上三角部の右に、下三角
部を圧縮した図、（ｆ）は列毎に、配列Ａに１次元的に
格納していく図である。 第２図は本発明を用いた時の行列×ベクトルのFORTRAN
コード例を示す図である。 第３図は行方向毎の格納法の説明図で、（ａ）は格納す
べき行列Ａの非ゼロ要素・ゼロ要素のパターン図、
（ｂ）は行列の非ゼロ要素格納の順序を示す図、（ｃ）
は同、対称行列の場合の非ゼロ要素格納の順序を示す
図、（ｄ）は配列Ａの先頭の数要素を示す図、（ｅ）は
各行の先頭の番地を配列するIAの内容を示す図、（ｆ）
は各要素の元の行列での列番号するJAの先頭の数要素を
示す図である。 第４図は行方向格納法を用いたときの行列×ベクトルの
FORTRANコード例を示す。 第５図は従来のベクトル計算機向きの格納法説明図で、
（ａ）は格納すべき行列の非ゼロ要素・ゼロ要素のパタ
ーン図、（ｂ）は同下三角部のゼロ要素・ゼロ要素のパ
ターン図、（ｃ）は同行列の上三角部のゼロ要素のパタ
ーン図、（ｃ）は同行列の上三角部のゼロ要素・ゼロ要
素のパターン図、（ｄ）は２次元配列に上三角部を詰め
た図、（ｅ）は同下三角部を詰めた図である。 第６図は第５図の格納法を用いたときの行列×ベクトル
のFORTRANコード例を示す。 第７図は線形ソートのアルゴリズムを示すフローチャー
ト図である。 第８図は行方向格納法から本発明の格納法へ変換するア
ルゴリズムを示すフローチャート図である。 第９図は有限要素法中における本発明の使用例のフロー
チャート図である。 第10図は実施例で用いた例題の領域分割図である。 第11図は各行列格納法を用いた行列×ベクトルの１要素
当たりの演算時間を示す図である。 第12図は各行列格納法を用いた、第９図の実施例に於け
るコストを示した図である。

Claims (1)

【特許請求の範囲】
1. 【請求項１】予めベクトルの番号を行毎の非ゼロ要素が
   単調減少になるよう付け替え、対角部分は別に格納し、
   非対角部は、上三角、下三角の順に左方に圧縮し列方向
   に一次元化して格納することを特徴とするベクトル計算
   機の記憶装置への大規模不規則疎行列格納法。