行列方程式計算装置および行列方程式計算方法
【課題】電磁界解析などの数値シミュレーションにおいて最終的に解くべき方程式として現れる大型疎行列方程式の計算を高速に行うことができ、しかもこの大型疎行列方程式の大型疎行列が対称行列である場合および非対称行列である場合の双方に適用することができる行列方程式計算装置および行列方程式計算方法を提供する。
【解決手段】大型疎行列方程式を反復解法により解く行列方程式計算装置１は、その大型疎行列方程式の大型疎行列の成分のうち、非ゼロの成分のみを１列ごとに格納する所定の個数のメモリを有するメモリ部と、反復解法の演算の少なくとも一部をデータフロー形式で実行する１つまたは複数の演算部とを有する。メモリ部から１行ごとの反復解法の演算に必要なデータを演算部に一度にロードするように構成する。
【解決手段】大型疎行列方程式を反復解法により解く行列方程式計算装置１は、その大型疎行列方程式の大型疎行列の成分のうち、非ゼロの成分のみを１列ごとに格納する所定の個数のメモリを有するメモリ部と、反復解法の演算の少なくとも一部をデータフロー形式で実行する１つまたは複数の演算部とを有する。メモリ部から１行ごとの反復解法の演算に必要なデータを演算部に一度にロードするように構成する。
【発明の詳細な説明】
【技術分野】
【０００１】
この発明は、行列方程式計算装置および行列方程式計算方法に関し、例えば、電磁界解析、流体解析、構造解析などの数値シミュレーションにおいて最終的に解くべき方程式として現れる大型（大規模）疎行列方程式の計算に適用して好適なものである。
【背景技術】
【０００２】
電磁界解析、流体解析、構造解析などにおける離散化された偏微分方程式を解くための科学技術計算は、多くの場合、最終的には、求めるべき未知数に関し行列方程式を解く問題に帰着する。また、これらの科学技術計算においては、多くの場合、解くべき行列方程式の行列は、次数が数千〜数千万と大規模で、かつ、その成分（要素）がほとんどゼロとなる疎行列、すなわち大型疎行列となる。実際、例えば、電磁界解析では、図１０に示すように、静電場・静磁場、準静的場、高周波の場などの電磁場の種類ごとに解くべき行列方程式は異なり、一方、偏微分方程式の解法には、大別して、差分法（ＦＤＭ）、有限要素法（ＦＥＭ）、境界要素法（ＢＥＭ）などの数値解析法があるため、これらの組合せにより非常に多くの計算方法が存在するものの、そのほとんどの計算方法は大型疎行列方程式を解く問題に帰着する。
【０００３】
これらの科学技術計算では、その計算規模が大きくなればなるほど全体の計算時間に対する行列方程式の計算時間の占める割合は大きくなり、しばしば９５％以上にもなる。近年のスーパーコンピュータの性能向上、および、パーソナルコンピュータ（ＰＣ）レベルにおける高性能化に助長され、科学技術計算に対する要求もますます強くなり、取り扱うべき問題も大規模化している。とりわけ、短期間での製品開発が要求される産業応用の場では、ハイパフォーマンスコンピューティング（ＨＰＣ）技術が非常に大きな解決課題となっている。このような背景から、科学技術計算の大部分を占める行列方程式の計算の高速化が、科学技術計算全体の最大の課題の１つとなっているといって過言でない。
【０００４】
一般に、行列方程式の解法は、原理的には未知数を１つずつ消去しながら解を求めていくＬＵ分解法に代表されるような直接解法と、初期値を設定してそこから真の解に向かって繰返し近似解を更新し収束値を求めていく反復解法とに大別される。このうち反復解法は、共役傾斜法（ＣＧ法）などのように線形連立方程式を等価なＮ次元２次形式（Ｎは行列の次数）の極小値を求める問題に置き換える方法と、そのような問題の置き換えは行わず行列を対角成分と非対角成分とに分離し直接対角成分のみから近似解の更新を行う単純反復法とに分類される。さらに、線形連立方程式を等価なＮ次元２次形式の極小値を求める問題に置き換える方法には、対称行列を適用対象とするＣＧ法をベースに、それを非対称行列に拡張した双共役傾斜法（ＢｉＣＧ法）や、反復計算を行う前に行列を前処理したり、各反復過程で最適な探索方向ベクトルを構成するなどして反復計算を安定化したり高速化したりする、ＩＣＣＧ法、ＣＧＳ法、ＢｉＣＧ−Ｓｔａｂ法、ＢｉＣＧ−Ｓｔａｂ２法などの派生的な方法が多数提案されている。
【０００５】
従来、これらの方法のうち大型疎行列方程式の解法としては、上記の反復解法が有効であると考えられており、電磁界解析、流体解析、構造解析などの各種の数値シミュレーションに現れる大型疎行列方程式の計算にしばしば用いられている。
【０００６】
ＣＧ法などの行列方程式の反復解法においては、解くべき行列形式の線形（１次）連立方程式の近似解を次の手順で求める。まず、Ｎを未知数の数あるいは行列の次数として、解くべき線形連立方程式を、それと等価なＮ次元２次形式の極小値を求める問題に置き換える。そして、初期値として適当なＮ次元解ベクトルを設定し、この初期値（ｋ＝０）を出発点にして、ＣＧ法などのそれぞれの反復解法の方法にしたがって次の近似解（ｋ＝１）を求め、さらに、この近似解（ｋ＝１）からその次の近似解（ｋ＝２）を求めていき、このｋ番目の近似解から（ｋ＋１）番目の近似解を求める手続きを繰返す。そして、この繰返し計算を（ｋ＋１）番目の近似解があらかじめ設定された精度以内になるまで行い、この精度以内になったところ、すなわち解が収束したところで上記の反復計算を終了し、この収束した値を数値計算としての最終的な近似解とする。
【０００７】
上記の行列方程式の反復解法の高速化技術としては、ＰＣクラスタ、あるいはさらに大規模なスーパーコンピュータなどの並列処理かつ高速なＣＰＵ上での走行が最も一般的に利用されているが、これらの計算機は比較的限られた研究機関でのみ利用が可能であったり、共同利用型であったりと、製品開発での使用などのだれもが手軽に利用できる一般の環境とはなっていない。
【０００８】
これに対し、近年のＦＰＧＡ、ＣＰＬＤなどの書換え可能なＬＳＩの出現およびその手軽な設計・利用環境、さらに、プリント回路基板設計・開発環境の充実化を背景に、科学技術計算のＨＰＣ技術の選択肢の一つとして、ターゲットをしぼり、計算方法に特化したアーキテクチャのハードウェアを構成することにより、安価に、かつ、身近なＰＣの環境で実現することができる専用計算機が新たに注目されつつある。
【０００９】
現在、線形行列計算を行うための専用計算機（以下、行列計算専用計算機ともいう）には、密行列をターゲットとしたものとして、ＬＵ分解法を高速に処理すべく複雑なＬＵ分解処理はホスト計算機に行わせ前進代入処理などの部分的な単純な処理を専用のハードウェアで行うタイプ、疎行列をターゲットとしたものとして、ハードウェア構成が比較的簡単な単純反復法の反復処理と収束判定処理とを直接ハードウェア化したタイプ、および、疎行列用のＬＵ分解法あるいはＣＧ法の演算のうち、ベクトルの内積演算あるいは行列とベクトルとの乗算の部分をハードウェア化したタイプなどがある（非特許文献１〜１３参照。）。
【００１０】
【非特許文献１】浅井秀樹，浅井光男，田中 衛，大規模スパース行列のＬＵ分解専用並列計算機，信学論，Vol.J69-D, No.7, pp.1044-1053 (1986)
【非特許文献２】田辺 昇，土肥康孝，連想スイッチによる疎行列計算機の構成，信学論，Vol.J70-D, No.12, pp.2393-2401 (1987)
【非特許文献３】田辺 昇，石坂健一，村越英機，泉谷昭二，土肥康孝，並列パイプライン式疎行列専用計算機ＲＡＭＰの構成，信学論，Vol.J71-D, No.10, pp.1939-1948 (1988)
【非特許文献４】清木 泰，福重俊幸，泰地真弘人，牧野淳一郎，小河正基，戒崎俊一，密行列専用計算機ＧＥＮＥＲＡＬ−１の開発，計算機アーキテクチャ，111-9, pp.65-72 (1995)
【非特許文献５】泰地真弘人，藤田 実，科学技術計算に適した準並列ＣＰＵの開発に関する研究，ＮＥＤＯ 平成１２年度提案公募事業成果報告会 予稿集，提案公募:99S29-019-1 (2000)
【非特許文献６】大野洋介，戒崎俊一，泰地真弘人，小長谷明彦，行列専用計算機ＭＡＣＥ（M Atrix Computation Engine）の実効性能
【非特許文献７】長谷川大，回路シミュレーション用行列演算専用計算機の開発，信学技報（TECHNICAL REPORT OF IEICE ），CA2002-116, pp.51-56 (2003)
【非特許文献８】瀧口 善貴, 松岡 俊佑, 川口 秀樹, ラプラス/ ポアソン方程式差分法ソルバー専用計算機の開発に関する研究，平成１７年電気・情報関係学会北海道支部連合大会 (2005), No.121.
【非特許文献９】G.R.Morris, V.K.Prasanna, R.D.Anderson, A Hybrid Approach for Mapping Conjugate Gradient onto an FPGA-Augmented Reconfigurable Supercomputer, Proc. of the 14th Annual IEEE Symposium on Field-Programmable Custom Computing Machines (FCCM'06), Napa, California (2006)
【非特許文献１０】R. Strzodka, D. Goeddeke, Pipelined Mixed Precision Algorithms on FPGAs for Fast and Accurate PDE Solvers from Low Precision Components, Extended technical report of IEEE Proc. on Field-Programmable Custom Computing Machines (FCCM'06), Napa, California (2006)
【非特許文献１１】O. Maslennikow, V. Lepekha, A. Sergyienko, FPGA Implementation of the Conjugate Gradient Method, LNCS 3911, pp.526-533 (2006)
【非特許文献１２】田原一平, 川口 秀樹, 松岡俊祐，線形計算専用計算機・共役傾斜法マシンの開発に関する研究，平成１９年電気・情報関係学会北海道支部連合大会 (2007), p.147.
【非特許文献１３】D. DuBois, A. DuBois, Sparse Matrix-Vector Multiplicationand Conjugate Gradient Algorithm on Hybrid Computing Platforms, LALP-07-041 (2007)
【発明の開示】
【発明が解決しようとする課題】
【００１１】
行列方程式計算用の専用計算機においては、ソフトウェアプログラムベースの計算機で行う処理と全く同じことをハードウェア化することは困難である。このため、従来は、全体の処理の中でハードウェア化に適し、かつ、ハードウェア化により有効に高速化が図れる処理を適切に選び出してそれらを専用計算機で実行させ、ホスト計算機と連動しながら全体として処理を高速化する方法が採られてきた。
【００１２】
非特許文献４では、密行列をターゲットとしガウスの消去法のための専用計算機を提案している。とりわけ、ガウスの消去法の中で演算量のほとんどを占める前進消去の内積演算の部分をハードウェア化する方式について検討し、計算機のシステム構成を提案している。しかしながら、メモリ容量の制約から、１０００次元以下の行列を対象としており、大規模な計算を行うにはさらなる検討が必要である。
【００１３】
非特許文献５では、同じく密行列をターゲットとしたＬＵ分解法の専用計算機について、ＬＵ分解法の計算ではデータ量（行列の次元）Ｎに対して、演算量がＮ2 、Ｎ1.5 とメモリアクセスがボトルネックになりづらい性質に着目して、ＬＵ分解法に現れる内積演算に特化し並列に高速処理することができる専用計算機の方式検討を行い、計算機のシステム構成を提案している。
【００１４】
また、非特許文献６では、同じく密行列をターゲットとしたＬＵ分解法の専用計算機について、ＬＵ分解法の演算のうち行列のＬＵ分解のための行列成分の四則演算処理のみをハードウェア化し、ピボット選択や前進代入・後退代入などの処理はホスト計算機で処理する方式について検討し、計算機のシステム構成を提案している。
【００１５】
このとき、上記の密行列を対象とした専用計算機による行列計算の高速化の試みにおいては、必然的に大容量のメモリが必要となり、また、より高速計算を実現するためにはハードウェアが大規模となるため、高々、数千次元程度の行列サイズを取り扱うのが限界であり、科学技術計算にしばしば現れるより大規模な計算への適用は困難であるという問題があった。
【００１６】
一方、非特許文献１〜３では、疎行列をターゲットとし疎行列用ＬＵ分解法のための専用計算機を提案している。疎行列性を利用し、行列の次元Ｎに比してはるかに少ない数の非ゼロ行列成分の分だけのローカルメモリを有する演算ユニットを用意し、１行分の加減算・乗除算を並列に行える構成とし、ＬＵ分解法の全ての処理を実行する方式について検討し、計算機のシステム構成を提案している。
【００１７】
また、非特許文献７では、同じく疎行列をターゲットとしたＬＵ分解法の専用計算機について、汎用の算術演算器を２つ用意し、かつ、メモリアクセスを２チャネルとした上でパイプライン処理が可能なハードウェアを小規模ながら実際のＰＣＩインターフェース機能を有するＦＰＧＡ評価キットに実装し、その動作確認まで行っている。
【００１８】
しかしながら、上記のＬＵ分解法をハードウェア化する試みにおいては、ピボット選択も含めた複雑な行列のＬＵ分解、そして、それに続く前進および後退代入の複数の処理が必要であり、ＣＧ法やＢｉＣＧ法などの反復解法に比べて演算量の多い計算を余儀なくされるという問題があった。
【００１９】
非特許文献８では、ラプラス・ポアソン方程式に特化した専用計算機を提案している。この専用計算機では、ラプラス・ポアソン方程式の差分式に現れる大型疎行列の１行ごとの非ゼロ成分の数が高々５であり、かつ、同じ行に現れる行列成分が隣り合うグリッドのアドレスとなることを利用し、この構造に特化した単純反復法のハードウェア化が提案されている。しかし、全体のシステムとしてシンプルなアーキテクチャで実現することができることが判明したものの、やはり、単純反復法は解の収束が悪く、行列解法のＨＰＣ化としては必ずしも適切なものとはなっていなかった。
【００２０】
非特許文献９では、大型疎行列をターゲットとしＣＧ法のための専用計算機を提案している。とりわけ、ＣＧ法の中で行列と探索方向のベクトルとの乗算演算の部分のみをハードウェア化し、他の処理はホスト計算機で実行させ、ホスト計算機と専用計算機とを連動させながらＣＧ法の計算を高速化する方式について検討し、計算機のシステム構成を提案している。
【００２１】
非特許文献１０では、大型疎行列をターゲットとしたＣＧ法の専用計算機について、ＦＰＧＡへの実装を想定した場合の精度とハードウェアサイズ・計算速度のトレードオフについて検討し、パイプライン処理を併用した適切な演算方法を提案している。
【００２２】
非特許文献１１では、大型疎行列をターゲットとしたＣＧ法の専用計算機について、ＣＧ法の演算に現れる除算の効率的な計算方法について検討し、演算部分について、実際にＦＰＧＡに実装しながらその有効性を示している。
【００２３】
非特許文献１２では、大型疎行列をターゲットとしたＣＧ法の専用計算機について、ＣＧ法の１回の反復計算のハードウェア化について検討し、小規模ながら実際にＦＰＧＡに実装しながらその実現可能性を示している。
【００２４】
非特許文献１３では、大型疎行列をターゲットとしたＣＧ法の専用計算機について、商用のリコンフィギュレーション可能な計算機を用いて行列方程式計算のＨＰＣ化について検討し、ＣＧ法に特化したシステムを提案している。
【００２５】
以上のように、従来の行列計算専用計算機は、必ずしも科学技術計算にしばしば現れる大型疎行列には適していないＬＵ分解法をハードウェア化するものや、あるいは大型疎行列に適したＣＧ法のハードウェア化である場合も、ＣＧ法の演算に除算などの複雑な処理が含まれるために、演算方法を部分的にハードウェア化し、ホスト計算機と連動して動作させる構成のものである。これらの行列計算専用計算機は、いずれも、完全な専用計算機という形態ではなく、ホスト計算機とのデータのやり取りを生じるなど動作速度、大規模計算化などの点で必ずしも十分なものではなかった。また、ＣＧ法は基本的に対称行列が対象であるが、科学技術計算では非対称行列もよく現れるため、対称行列に特化した専用計算機では適用範囲が限られてしまうという問題もあった。
【００２６】
そこで、この発明が解決しようとする課題は、電磁界解析、流体解析、構造解析などの各種の数値シミュレーションにおいて最終的に解くべき方程式として現れる大型疎行列方程式の計算を高速に行うことができ、しかもこの大型疎行列方程式の大型疎行列が対称行列である場合および非対称行列である場合の双方に適用することができる行列方程式計算装置および行列方程式計算方法を提供することである。
【課題を解決するための手段】
【００２７】
上記課題を解決するために、第１の発明は、
対象となる大型疎行列方程式を反復解法により解く行列方程式計算装置であって、
上記大型疎行列方程式の大型疎行列の成分のうち、非ゼロの成分のみを１列ごとに格納する所定の個数のメモリを有するメモリ部と、
上記反復解法の演算の少なくとも一部をデータフロー形式で実行する１つまたは複数の演算部とを有し、
上記メモリ部から１行ごとの反復解法の演算に必要なデータを上記演算部に一度にロードするように構成されていることを特徴とするものである。
この行列方程式計算装置は、典型的には、パーソナルコンピュータなどからなるホスト計算機と接続され、このホスト計算機から必要なデータを取得するが、これに限定されるものではなく、必要に応じて、ホスト計算機の機能を持たせてもよい。
【００２８】
第２の発明は、
対象となる大型疎行列方程式を反復解法により解く行列方程式計算方法であって、
上記大型疎行列方程式の大型疎行列の成分のうち、非ゼロの成分のみを１列ごとにメモリ部の所定の個数のメモリに格納するステップと、
上記メモリ部から１行ごとの反復解法の演算に必要なデータを１つまたは複数の演算部に一度にロードするステップと、
上記反復解法の演算の少なくとも一部を上記演算部でデータフロー形式で実行するステップとを有することを特徴とするものである。
【００２９】
第１および第２の発明において、大型疎行列とは、電磁界解析、流体解析、構造解析などの数値シミュレーションにおいて最終的に解くべき行列方程式に現れる数千〜数百万あるいはそれ以上の次元であるがその成分のほとんどはゼロであり、非ゼロ成分は、１行あたり数個〜十数個であるような行列を意味する。大型疎行列は対称行列であっても非対象行列であってもよい。反復解法には、例えば、ＣＧ法、ＢｉＣＧ法、ＩＣＣＧ法、ＢｉＣＧ−Ｓｔａｂ法、ＣＧＳ法、ＢｉＣＧ−Ｓｔａｂ２法などを用いることができる。
【００３０】
第１および第２の発明においては、大型疎行列の成分のうち、非ゼロの成分のみを１列ごとに格納する所定の個数（複数）（必要に応じて選ばれるが、例えば８あるいは１６）のメモリを用いるが、これは次のような理由による。図１は大型疎行列方程式の１例を示す。図１に示すように、行列Ａおよび非同次ベクトルｂからなる大型疎行列方程式
【数１】
から未知数ベクトルｘを求める問題において、これをＣＧ法の計算方法
【数２】
を用いて解く場合を考える。（２−２）式の反復部分の計算には、
【数３】
なる行列とベクトルとの乗算がある。このベクトルｑのｉ番目の成分ｑi の計算には行列Ａのｉ行目の１行分の成分の値が必要であるが、行列の成分の値を１つのメモリに格納しては、必要な値を１度にこれらのメモリから取得することはできない。このため、第１および第２の発明では、メモリ部が、行列の非ゼロ成分のみを１列ごとに格納する複数のメモリを有する構成とし、１行ごとの計算に必要なデータをメモリ部から演算部に並列に一度にロードすることとしている。また、後述の複数の演算部による並列処理の場合は、これらの列ごとにそれぞれ設けられたメモリをさらに並列数分設け、それぞれの並列処理が分担する行の成分の値を格納する。
【００３１】
好適には、行列方程式計算装置は、同じハードウェア回路の複数の演算部を有する。そして、これらの演算部の数だけ、行列の行数Ｎを分割し、それぞれの演算部が連動動作しながら並列処理を行う。こうすることで、行列方程式の計算を演算部の数だけスケーラブルに並列処理により高速化することができる。
好適には、反復解法の演算を可能な限りデータフロー形式で実行する。
【００３２】
さらに、好適には、対称行列を対象としたＣＧ法専用計算機の場合は、偶数かつ複数の同じハードウェア回路の演算部を有する。そして、これらの演算部の２つをペアに連動させることにより、１つの非対称行列を対象としたＢｉＣＧ法専用計算機としても機能させる。また、この演算部のペアの数だけ行列の行数を分割し、それぞれの演算部が連動動作しながら並列処理を行うようにしてもよい。こうすることで、１つの専用計算機を対称行列および非対称行列両方に適用することができ、対称行列モードの場合は非対称行列モードの場合の２倍の速度で動作する計算装置を実現することができる。
【発明の効果】
【００３３】
この発明によれば、行列方程式計算装置にＯＳなどのソフトウェア的な制約なしに、ハードウェア的に許される限り大容量のメモリを搭載可能となるため、大規模な計算を実行することができる。また、演算部におけるＣＧ法などの反復解法のデータフローアーキテクチャ回路構成とメモリ部における行列の１列の非ゼロ成分ごとのメモリ構成とにより、スループットとしても、例えば、ＣＧ法やＢｉＣＧ法では１回の反復計算を３Ｎ＋２クロック（Ｎは行列の次数）で、ＢｉＣＧ−Ｓｔａｂ法やＣＧＳ法では５Ｎクロックで処理することができ、計算方法の限界まで動作時間を最小化した高速な動作を実現することができる。
また、複数の演算部を有する構成とし並列動作させることにより、スケーラブルな並列処理で高速化を図ることができる。
さらに、ＣＧ法の演算部を２つ連動させることでＢｉＣＧ法の演算部を実現することができるため、用途に適した効率的な計算を行うことができる。
さらに、複数の反復解法の演算部を用意し、これらの演算部に同じ行列計算を同時に行わせ、最も早く計算が完了した演算部の解をホスト計算機がアップロードすることにより、自動的に問題に最適な反復解法の方法を選択して行列方程式の計算を実行することができる。
【発明を実施するための最良の形態】
【００３４】
以下、この発明の実施の形態について図面を参照しながら説明する。
まず、この発明の第１の実施の形態による行列方程式計算装置について説明する。
この行列方程式計算装置はＣＧ法専用計算機により構成される。この行列方程式計算装置は、（１）式の行列方程式を（２−１）式、（２−２）式および（３）式に示すＣＧ法により計算するものである。
【００３５】
図２はこの行列方程式計算装置１の全体構成、図３はこの行列方程式計算装置１の各部の構成の詳細を示す。図２および図３に示すように、この行列方程式計算装置１は、大型疎行列Ａの全成分の値を格納するメモリモジュール１１と、（２−２）式のＣＧ法の１回分の反復計算をデータフロー形式で３Ｎクロックで実行する演算モジュール１２と、メモリモジュール１１と演算モジュール１２との間で、（３）式の大型疎行列Ａとベクトルｐk との乗算をクロック遅延のない組合せ回路で実行し、この結果の値Ａｐk を演算モジュール１２に提供する行列−ベクトル乗算回路１３と、演算モジュール１２で得られたｋ＋１番目の反復解が収束しており最終的な解となっているかどうかを判定する収束判定部１４とを有する。メモリモジュール１１、行列−ベクトル乗算回路１３、演算モジュール１２および収束判定部１４における一連のデータの流れの動作はマスターコントローラ１５により制御される。
【００３６】
メモリモジュール１１は、大型疎行列Ａを非ゼロ成分のみに圧縮したＮ×ｍのサイズの行列Ａ’の成分を列ごとに格納するｍ個のメモリＭ１−１〜Ｍ１−ｍを有する。この大型疎行列Ａの圧縮は具体的には次のように行う。例えば、各行ごとの非ゼロ成分の数のうち全Ｎ行の中で最大のものをｍとしてＮ×ｍのサイズの行列を用意し、これに各行ごとに非ゼロ成分を左詰めで格納する。同時に、これらのメモリＭ１−１〜Ｍ１−ｍと全く同じサイズのメモリをもう１つずつ用意する（メモリＭ２−１〜Ｍ２−ｍ）。これらのメモリＭ２−１〜Ｍ２−ｍには、対応するメモリＭ１−１〜Ｍ１−ｍの非ゼロ成分がもともとの大型疎行列Ａのどの列に位置していたかを示すインデックスを格納する。これにより、通常、大型疎行列ではＮが数千〜数千万に対しｍは高々十数個であるため、大型疎行列Ａの情報を失うことなく大幅にメモリ容量を節約することができる。メモリモジュール１１をこのように構成することで、メモリ容量を大幅に節約することができ、かつ、一度に１行分の行列の非ゼロ成分全てに並列アクセスすることができるためスループットとしての計算性能に大きく影響する行列の成分が格納されたメモリＭ１−１〜Ｍ１−ｍへのアクセスは、最小限の１回にまで減らすことができる。
【００３７】
メモリモジュール１１はさらに、ＣＧ法の反復計算実行時に一時的に値を格納しておくために、同じベクトルｐk の値をｍ個別々に格納するためのメモリＭ３−０〜Ｍ３−ｍ、ベクトルｐk の更新値ｐk+1 を格納するためのメモリＭ３−０’、ベクトルｘk を格納するためのメモリＭ４、ベクトルｘk の更新値ｘk+1 を格納するためのメモリＭ４’、ベクトルｒk を格納するためのメモリＭ５、ベクトルｒk の更新値ｒk+1 を格納するためのメモリＭ５’、ベクトルｑk を格納するためのメモリＭ６、ベクトルｒk のノルムの２乗‖ｒk ‖2 ＝（ｒk ，ｒk ）を格納するためのメモリＭ７、ノルムの２乗‖ｒk ‖2 の更新値‖ｒk+1 ‖2 ＝（ｒk+1 ，ｒk+1 ）を格納するためのメモリＭ７’、行列の次数Ｎ、行ごとの非ゼロ成分の最大数ｍ、非同次項ベクトルｂのノルムの２乗‖ｂ‖2 ＝（ｂ，ｂ）、相対残差の収束判定基準値ε、最大反復回数Ｋｍ、各反復ステップの残差Ｒｓなどの反復の収束条件を格納するためのメモリＭ８を有する。
【００３８】
上述のメモリＭ１−１〜Ｍ１−ｍ、Ｍ２−１〜Ｍ２−ｍ、Ｍ３−０〜Ｍ３−ｍ、Ｍ３−０’、Ｍ４、Ｍ４’、Ｍ５、Ｍ５’、Ｍ６、Ｍ７、Ｍ７’、Ｍ８としては、それぞれ独立した個別のメモリを用いてもよいし、それらのうちの１つまたは２つ以上のメモリとして、独立した１つのメモリのメモリ領域を所定個数に分割したものを用い、これらのメモリ領域をメモリＭ１−１〜Ｍ１−ｍ、Ｍ２−１〜Ｍ２−ｍ、Ｍ３−０〜Ｍ３−ｍ、Ｍ３−０’、Ｍ４、Ｍ４’、Ｍ５、Ｍ５’、Ｍ６、Ｍ７、Ｍ７’、Ｍ８のうちのいずれかに用いてもよい。これらのメモリＭ１−１〜Ｍ１−ｍ、Ｍ２−１〜Ｍ２−ｍ、Ｍ３−０〜Ｍ３−ｍ、Ｍ３−０’、Ｍ４、Ｍ４’、Ｍ５、Ｍ５’、Ｍ６、Ｍ７、Ｍ７’、Ｍ８としては、例えば、ＳＲＡＭまたはＤＲＡＭを用いることができる。
【００３９】
行列−ベクトル乗算回路１３は、上述のように列ごとにメモリＭ１−１〜Ｍ１−ｍに分割格納された大型疎行列Ａの成分を、メモリＭ２−１〜Ｍ２−ｍを参照しつつ１行ごとに１行の全ての成分に並列に同時アクセスし、また、ベクトルｐk を格納するメモリＭ３−１〜Ｍ３−ｍからそれぞれ（３）式の計算に必要な値のアドレスに同時アクセスし、データフロー形式で（３）式の計算を実行し、ベクトルｑk を１クロックで１成分ずつ計算し、その結果を演算モジュール１２に引き渡す。
【００４０】
演算モジュール１２は、行列−ベクトル乗算回路１３からのベクトルｑk 、メモリＭ３−０のベクトルｐk 、メモリＭ４のベクトルｘk 、メモリＭ５のベクトルｒk 、メモリＭ６のベクトルｑk およびメモリＭ７の（ｒk ，ｒk ）を用いて、（２−２）式（ｉ）中のベクトルｐk とベクトルｑk （＝Ａｐk ）との内積（ｐk ，ｑk ）の計算を実行するための演算回路３１、（２−２）式（ｉ）中の除算（ｒk ，ｒk ）／（ｐk ，ｑk ）の計算を実行するための演算回路３２、（２−２）式（ｉｉ）のベクトルｘk+1 の計算を実行するための演算回路３３、（２−２）式（ｉｉｉ）のベクトルｒk+1 の計算を実行するための演算回路３４、（２−２）式（ｖ）中のベクトルｒk+1 の２乗（ｒk+1 ，ｒk+1 ）の計算を実行するための演算回路３５、（２−２）式（ｉｖ）中の除算（ｒk+1 ，ｒk+1 ）／（ｒk ，ｒk ）の計算を実行するための演算回路３６、（２−２）式（ｖ）のベクトルｐk+1 の計算を実行するための演算回路３７を有する。
【００４１】
（２−２）式からわかるように、演算回路３１、３３、３４、３５、３７の計算はそれぞれＮクロックで、演算回路３２、３６の計算は１クロックで実行することができるが、演算モジュール１２は全体として可能な限りこれらの計算を並列に実行することができるように構成されている。具体的には、演算回路３１、３２のＮ＋１クロックの計算が完了した後は、演算回路３３、３４、３５の計算は、演算回路３４、３５のデータフロー性を利用するとまとめてＮクロックで実行することができるため、残りの演算回路３６、３７のＮ＋１クロックの計算を合わせて全体で３Ｎ＋２クロックで実行することができるように構成されている。
【００４２】
収束判別部１４は、ＣＧ法の演算を実行する行列−ベクトル乗算回路１３から演算モジュール１２での１回分の反復計算で得られたｋ＋１回目の残差ベクトルｒk+1 から相対誤差を計算し、あらかじめ設定されたメモリＭ８内の相対残差の収束判定基準値ε内に収まり反復計算が収束しているかどうかを判別する。
【００４３】
マスターコントローラ１５は、計算が収束していればｘk+1 を最終的な近似解として計算を終了し、ホスト計算機（図示せず）に正常に収束値を得た旨を通知する。計算が収束していなければ、あらかじめ設定されたメモリＭ８内の最大反復回数Ｋｍを超えた場合は、ホスト計算機に反復計算が収束しなかった旨を通知し、最大反復回数Ｋｍを超えていない場合は、反復計算で得られたベクトルｐk+1 、ｘk+1 、ｒk+1 および（ｒk+1 ，ｒk+1 ）を次の反復計算に送り、再度、行列−ベクトル乗算回路１３から演算モジュール１２での反復計算の実行を指示する。
【００４４】
図４は行列方程式計算装置１の使用形態を示す。行列方程式計算装置１はＰＣＩインターフェース４１によりホスト計算機４２と接続されている。ＣＧ法の演算を実行する上で必要な値は、圧縮格納された大型疎行列Ａ、反復計算の初期値となるベクトルｐ0 （＝ｒ0 ）、ｘ0 、非同次項の（ｂ，ｂ）、収束判定基準ε、最大反復回数Ｋｍである。これらの値はホスト計算機４２に格納され、計算開始指示の前に行列方程式計算装置１にダウンロードされてメモリモジュール１１のメモリＭ３、Ｍ４、Ｍ５、Ｍ８に格納されるようになっている。ホスト計算機４２は、これらの値を行列方程式計算装置１にダウンロードした後、行列方程式計算装置１に対し計算開始を指示する。行列方程式計算装置１は、（２−２）式の反復計算および収束判定を繰返し、最大反復回数Ｋｍ以内で解が収束判定基準ε以下になれば、メモリＭ４’のベクトルｘk+1 を近似解としてホスト計算機４２に引渡し、最大反復回数Ｋｍまで計算しても収束しなければ、計算は未収束として、これをホスト計算機４２に通知する。
【００４５】
図５Ａ〜Ｅは行列方程式計算装置１の実装例を示す。図５ＡおよびＢに示すように、この例では、１つの小プリント基板５１上に、メモリモジュール１１の２列分の圧縮行列の成分を格納するメモリ５２と、この２列の行列と対応する探索ベクトルｐk の成分との乗算を行う行列−ベクトル乗算回路１３とが搭載される。図５ＣおよびＤに示すように、この小プリント基板５１が合計１６個、サブプリント基板上５３に接続される。図５Ｅに示すように、メインプリント基板上５４に、演算モジュール１２、収束判定部１４およびマスターコントローラ１５が搭載され、これにサブプリント基板５３が接続される。こうして、ＣＧ法専用計算機としての行列方程式計算装置１が構成される。これらの搭載部品のほとんどは例えばＦＰＧＡまたはＡＳＩＣにより作製することができる。このように構成することにより、行列方程式計算装置１をある決められた列数ｍ用に構成した場合でも、あらかじめ十分なソケットを用意しておけば、後に小プリント基板５１を追加接続することにより、行列サイズの変更に対してフレキシビリティーを持たすことができる。
【００４６】
この第１の実施の形態によれば、次のような利点を得ることができる。すなわち、行列方程式計算装置１をＣＧ法専用計算機により構成し、メモリモジュール１１と演算モジュール１２とを分離した構成としていることにより、数値シミュレーションに用いる大型疎行列に応じてメモリモジュール１１を適切な構成のものに交換することができ、比較的複雑な構成を有する演算モジュール１３は変更しないで、大型疎行列の変更に柔軟に対応することができ、フレキシビリティーが高い。また、メモリモジュール１１は、行列の１列ごとにｍ個のメモリＭ１−１〜Ｍ１−ｍを設け、これらのメモリＭ１−１〜Ｍ１−ｍから（３）式の１行分の計算に必要なデータを１クロックで行列−ベクトル乗算回路１３にロードし、（３）式の１行分の計算を並列的に実行するようにしていることにより、メモリモジュール１１へのアクセス回数を最小限にすることができ、それによって高速動作化を図ることができ、行列方程式の計算を高速かつ短時間で行うことができる。また、この行列方程式計算装置１は実装や計算の大規模化が容易である。以上により、極めて実用性が高い高性能の行列方程式計算装置１を実現することができる。
この行列方程式計算装置１は、電磁界解析、流体解析、構造解析などの各種の数値シミュレーションで最終的に解くべき方程式として現れる大型疎行列方程式を解くのに適用して好適なものである。
【００４７】
次に、この発明の第２の実施の形態による行列方程式計算装置について説明する。
この第２の実施の形態においては、メモリモジュール１１へのアクセス回数の最小限化に加え、さらに、演算モジュール１２における３Ｎ＋２クロック分の計算をより高速化すべく、演算モジュール１２を複数用意し、３Ｎ＋２クロック分の計算をこれらの複数の演算モジュール１２で並列に実行する方式が用いられる。演算モジュール１２における３Ｎ＋２クロックの計算はもともと全て並列処理してもよく、この部分の計算を複数の演算モジュール１２で行うことにより、１つの演算モジュール１２を用いた場合に比べ、演算モジュール１２の台数分だけ、高速化を図ることができる。上記以外のことは第１の実施の形態と同様である。
この行列方程式計算装置１の実装は、図６に示すように、図５Ｅに示すものと同様なメインプリント基板５４上にサブプリント基板５３を複数搭載することにより簡単に行うことができる。
【００４８】
この第２の実施の形態によれば、第１の実施の形態と同様な利点に加えて、動作速度のより一層の向上を図ることができ、行列方程式の計算をより高速に行うことができるという利点を得ることができる。
【００４９】
次に、この発明の第３の実施の形態による行列方程式計算装置について説明する。
この第３の実施の形態においては、複数の演算モジュール３の並列処理による高速化に加え、さらに、ＣＧ法は大型疎行列が対称行列である場合のみに適用が限定されているのに対し、非対称行列である場合にも適用することができる機能を装備すべく、２つの演算モジュール１２を連動動作させ、非対称行列にも適用することができるＢｉＣＧ法マシンとして動作させる方式が用いられる。ＢｉＣＧ法はもともとＣＧ法に転置行列に対するＣＧ法の処理を追加したものであり、もう１つの演算モジュール１２に転置行列を処理させることにより、２つの演算モジュール１２をペアで連動動作させることによりＢｉＣＧ法の動作を実現することができる。上記以外のことは第１の実施の形態と同様である。
この行列方程式計算装置１の実装は、図６に示すように、図５Ｅに示すものと同様なメインプリント基板５４上にサブプリント基板５３を複数搭載することにより簡単に行うことができる。
【００５０】
次に、ハードウェア記述言語であるＶＨＤＬ（VHSIC Hardware Description Language)により行列方程式計算装置１による大型疎行列計算の論理シミュレーションを行った結果について説明する。
図７はシミュレーションの例として使用したリニアモータの界磁部のモデルを、図８はその数値モデルを示す。このモデルは、３つの励磁コイル６１〜６３で励磁された静磁場が、鉄（μr ＝５０００）からなる磁場遮蔽材料６４によって外部にもれないよう遮蔽される現象に関するもので、系は軸対称性を持っており、図８の数値シミュレーションの解析領域は軸対称性を考慮して２次元の５０×１００グリッドサイズのものである。この磁場を記述する支配方程式は、与えられた電流密度ベクトルＪ、透磁率μに対し、ベクトルポテンシャルＡを未知数とし、
【数４】
となる。この方程式を半径方向、軸方向ともにΔｌのサイズの一様なグリッドで分割した上で軸対称性を考慮して差分法に基づき差分式で表すと、半径方向にi 番目、軸方向にj
番目のグリッドでは、
【数５】
となる。ただし、μi,j 、Ａi,j 、Ｊi,j はそれぞれ、半径方向にi 番目、軸方向にj 番目のグリッドでの透磁率μ、ベクトルポテンシャルＡのθ成分、電流密度ベクトルＪのθ成分である。このとき、未知数Ａi,j の数は、全てのグリッド数の５０×１００＝５０００であるので、差分法により数値解析する際に現れる大型疎行列は５０００×５０００のサイズとなる。この大型疎行列計算を、非対称行列用のＢｉＣＧ法に基づいて構成した行列方程式計算装置１により実行する場合のＶＨＤＬ論理回路シミュレーションにより行った。その解から得られた静磁場分布の結果を図９に示す。実際の計算では、（５）式からわかるように、１行に現れる未知数Ａi,j の係数は高々５個であるので、この疎行列性を利用し、行列は５０００×５のサイズに圧縮したものを用いている。図９に示す結果は、同じ数値モデルに基づいてＣ言語ソフトウェアシミュレーションを行った結果と一致しており、この行列方程式計算装置１の妥当性が保証される。このＶＨＤＬ論理回路シミュレーションでの動作を参考に、この行列方程式計算装置１がアクセス速度２６６ＭＨｚのメモリと同程度の速度で動作した場合の性能を見積ると、２．６ＧＨｚのＣＰＵを有するパソコン上でＣ言語を用いて計算した時の約５０〜６０倍の性能に匹敵する。
【００５１】
ここでは、簡単のため、軸対称２次元となるモデルを用いたが、基本的には３次元の場合も全く同様である。すなわち、例えば同じ図７の３次元モデルでは、グリッドサイズが１００×１００×１００であるので、最終的に現れる大型疎行列は、１００００００×１００００００で、（４）式に対応する３次元の差分式では、係数は１６個となるため、最終的には、１００００００×１６のサイズに圧縮された行列を解くことに帰着する。このため、この行列方程式計算装置１では、軸対称２次元の場合が、Ｎ＝５０００、ｍ＝５であったものが、３次元では、単に、Ｎ＝１００００００、ｍ＝１６となるだけで、動作自体は全く同様である。
この第３の実施の形態によれば、第１の実施の形態と同様な利点に加えて、大型疎行列Ａが非対称行列である場合にも行列方程式計算装置１を適用することができるという利点を得ることができる。
【００５２】
以上、この発明の実施の形態について具体的に説明したが、この発明は、上述の実施の形態に限定されるものではなく、この発明の技術的思想に基づく各種の変形が可能である。
例えば、上述の実施の形態において挙げた数値、回路構成、配置などはあくまでも例に過ぎず、必要に応じて、これらと異なる数値、回路構成、配置などを用いてもよい。また、行列方程式の解法には、ＣＧ法やＢｉＣＧ法ではなく、ＩＣＣＧ法、ＢｉＣＧ−Ｓｔａｂ法、ＣＧＳ法、ＢｉＣＧ−Ｓｔａｂ２法などの他の反復解法を採用し、演算モジュール１２もこれらの行列方程式の解法用のもので構成してもよい。
【図面の簡単な説明】
【００５３】
【図１】大型疎行列方程式の１例を示す略線図である。
【図２】この発明の第１の実施の形態による行列方程式計算装置の全体構成を示す略線図である。
【図３】この発明の第１の実施の形態による行列方程式計算装置の各部の構成の詳細を示す略線図である。
【図４】この発明の第１の実施の形態による行列方程式計算装置をホスト計算機と接続した状態を示す略線図である。
【図５】この発明の第１の実施の形態による行列方程式計算装置の実装例を示す略線図である。
【図６】この発明の第２の実施の形態による行列方程式計算装置の実装例を示す略線図である。
【図７】この発明の第３の実施の形態による行列方程式計算装置を用いて行った数値シミュレーションに用いた数値モデルを示す略線図である。
【図８】この発明の第３の実施の形態による行列方程式計算装置を用いて行った数値シミュレーションに用いた数値モデルを示す略線図である。
【図９】この発明の第３の実施形態による行列方程式計算装置を用いて行った数値シミュレーションの結果を示す略線図である。
【図１０】電磁界解析を例に、多くの科学技術計算は最終的に大型疎行列方程式の計算に帰着することを説明するための略線図である。
【符号の説明】
【００５４】
１…行列方程式計算装置、１１…メモリモジュール、１２…演算モジュール、１３…行列−ベクトル乗算回路、１４…収束判定部、１５…マスターコントローラ、４２…ホスト計算機、５１…プリント基板、５３…サブプリント基板、５４…メインプリント基板、Ｍ１−１〜Ｍ１−ｍ、Ｍ２−１〜Ｍ２−ｍ、Ｍ３−０〜Ｍ３−ｍ、Ｍ３−０’、Ｍ４、Ｍ４’、Ｍ５、Ｍ５’、Ｍ６、Ｍ７、Ｍ７’、Ｍ８…メモリ
【技術分野】
【０００１】
この発明は、行列方程式計算装置および行列方程式計算方法に関し、例えば、電磁界解析、流体解析、構造解析などの数値シミュレーションにおいて最終的に解くべき方程式として現れる大型（大規模）疎行列方程式の計算に適用して好適なものである。
【背景技術】
【０００２】
電磁界解析、流体解析、構造解析などにおける離散化された偏微分方程式を解くための科学技術計算は、多くの場合、最終的には、求めるべき未知数に関し行列方程式を解く問題に帰着する。また、これらの科学技術計算においては、多くの場合、解くべき行列方程式の行列は、次数が数千〜数千万と大規模で、かつ、その成分（要素）がほとんどゼロとなる疎行列、すなわち大型疎行列となる。実際、例えば、電磁界解析では、図１０に示すように、静電場・静磁場、準静的場、高周波の場などの電磁場の種類ごとに解くべき行列方程式は異なり、一方、偏微分方程式の解法には、大別して、差分法（ＦＤＭ）、有限要素法（ＦＥＭ）、境界要素法（ＢＥＭ）などの数値解析法があるため、これらの組合せにより非常に多くの計算方法が存在するものの、そのほとんどの計算方法は大型疎行列方程式を解く問題に帰着する。
【０００３】
これらの科学技術計算では、その計算規模が大きくなればなるほど全体の計算時間に対する行列方程式の計算時間の占める割合は大きくなり、しばしば９５％以上にもなる。近年のスーパーコンピュータの性能向上、および、パーソナルコンピュータ（ＰＣ）レベルにおける高性能化に助長され、科学技術計算に対する要求もますます強くなり、取り扱うべき問題も大規模化している。とりわけ、短期間での製品開発が要求される産業応用の場では、ハイパフォーマンスコンピューティング（ＨＰＣ）技術が非常に大きな解決課題となっている。このような背景から、科学技術計算の大部分を占める行列方程式の計算の高速化が、科学技術計算全体の最大の課題の１つとなっているといって過言でない。
【０００４】
一般に、行列方程式の解法は、原理的には未知数を１つずつ消去しながら解を求めていくＬＵ分解法に代表されるような直接解法と、初期値を設定してそこから真の解に向かって繰返し近似解を更新し収束値を求めていく反復解法とに大別される。このうち反復解法は、共役傾斜法（ＣＧ法）などのように線形連立方程式を等価なＮ次元２次形式（Ｎは行列の次数）の極小値を求める問題に置き換える方法と、そのような問題の置き換えは行わず行列を対角成分と非対角成分とに分離し直接対角成分のみから近似解の更新を行う単純反復法とに分類される。さらに、線形連立方程式を等価なＮ次元２次形式の極小値を求める問題に置き換える方法には、対称行列を適用対象とするＣＧ法をベースに、それを非対称行列に拡張した双共役傾斜法（ＢｉＣＧ法）や、反復計算を行う前に行列を前処理したり、各反復過程で最適な探索方向ベクトルを構成するなどして反復計算を安定化したり高速化したりする、ＩＣＣＧ法、ＣＧＳ法、ＢｉＣＧ−Ｓｔａｂ法、ＢｉＣＧ−Ｓｔａｂ２法などの派生的な方法が多数提案されている。
【０００５】
従来、これらの方法のうち大型疎行列方程式の解法としては、上記の反復解法が有効であると考えられており、電磁界解析、流体解析、構造解析などの各種の数値シミュレーションに現れる大型疎行列方程式の計算にしばしば用いられている。
【０００６】
ＣＧ法などの行列方程式の反復解法においては、解くべき行列形式の線形（１次）連立方程式の近似解を次の手順で求める。まず、Ｎを未知数の数あるいは行列の次数として、解くべき線形連立方程式を、それと等価なＮ次元２次形式の極小値を求める問題に置き換える。そして、初期値として適当なＮ次元解ベクトルを設定し、この初期値（ｋ＝０）を出発点にして、ＣＧ法などのそれぞれの反復解法の方法にしたがって次の近似解（ｋ＝１）を求め、さらに、この近似解（ｋ＝１）からその次の近似解（ｋ＝２）を求めていき、このｋ番目の近似解から（ｋ＋１）番目の近似解を求める手続きを繰返す。そして、この繰返し計算を（ｋ＋１）番目の近似解があらかじめ設定された精度以内になるまで行い、この精度以内になったところ、すなわち解が収束したところで上記の反復計算を終了し、この収束した値を数値計算としての最終的な近似解とする。
【０００７】
上記の行列方程式の反復解法の高速化技術としては、ＰＣクラスタ、あるいはさらに大規模なスーパーコンピュータなどの並列処理かつ高速なＣＰＵ上での走行が最も一般的に利用されているが、これらの計算機は比較的限られた研究機関でのみ利用が可能であったり、共同利用型であったりと、製品開発での使用などのだれもが手軽に利用できる一般の環境とはなっていない。
【０００８】
これに対し、近年のＦＰＧＡ、ＣＰＬＤなどの書換え可能なＬＳＩの出現およびその手軽な設計・利用環境、さらに、プリント回路基板設計・開発環境の充実化を背景に、科学技術計算のＨＰＣ技術の選択肢の一つとして、ターゲットをしぼり、計算方法に特化したアーキテクチャのハードウェアを構成することにより、安価に、かつ、身近なＰＣの環境で実現することができる専用計算機が新たに注目されつつある。
【０００９】
現在、線形行列計算を行うための専用計算機（以下、行列計算専用計算機ともいう）には、密行列をターゲットとしたものとして、ＬＵ分解法を高速に処理すべく複雑なＬＵ分解処理はホスト計算機に行わせ前進代入処理などの部分的な単純な処理を専用のハードウェアで行うタイプ、疎行列をターゲットとしたものとして、ハードウェア構成が比較的簡単な単純反復法の反復処理と収束判定処理とを直接ハードウェア化したタイプ、および、疎行列用のＬＵ分解法あるいはＣＧ法の演算のうち、ベクトルの内積演算あるいは行列とベクトルとの乗算の部分をハードウェア化したタイプなどがある（非特許文献１〜１３参照。）。
【００１０】
【非特許文献１】浅井秀樹，浅井光男，田中 衛，大規模スパース行列のＬＵ分解専用並列計算機，信学論，Vol.J69-D, No.7, pp.1044-1053 (1986)
【非特許文献２】田辺 昇，土肥康孝，連想スイッチによる疎行列計算機の構成，信学論，Vol.J70-D, No.12, pp.2393-2401 (1987)
【非特許文献３】田辺 昇，石坂健一，村越英機，泉谷昭二，土肥康孝，並列パイプライン式疎行列専用計算機ＲＡＭＰの構成，信学論，Vol.J71-D, No.10, pp.1939-1948 (1988)
【非特許文献４】清木 泰，福重俊幸，泰地真弘人，牧野淳一郎，小河正基，戒崎俊一，密行列専用計算機ＧＥＮＥＲＡＬ−１の開発，計算機アーキテクチャ，111-9, pp.65-72 (1995)
【非特許文献５】泰地真弘人，藤田 実，科学技術計算に適した準並列ＣＰＵの開発に関する研究，ＮＥＤＯ 平成１２年度提案公募事業成果報告会 予稿集，提案公募:99S29-019-1 (2000)
【非特許文献６】大野洋介，戒崎俊一，泰地真弘人，小長谷明彦，行列専用計算機ＭＡＣＥ（M Atrix Computation Engine）の実効性能
【非特許文献７】長谷川大，回路シミュレーション用行列演算専用計算機の開発，信学技報（TECHNICAL REPORT OF IEICE ），CA2002-116, pp.51-56 (2003)
【非特許文献８】瀧口 善貴, 松岡 俊佑, 川口 秀樹, ラプラス/ ポアソン方程式差分法ソルバー専用計算機の開発に関する研究，平成１７年電気・情報関係学会北海道支部連合大会 (2005), No.121.
【非特許文献９】G.R.Morris, V.K.Prasanna, R.D.Anderson, A Hybrid Approach for Mapping Conjugate Gradient onto an FPGA-Augmented Reconfigurable Supercomputer, Proc. of the 14th Annual IEEE Symposium on Field-Programmable Custom Computing Machines (FCCM'06), Napa, California (2006)
【非特許文献１０】R. Strzodka, D. Goeddeke, Pipelined Mixed Precision Algorithms on FPGAs for Fast and Accurate PDE Solvers from Low Precision Components, Extended technical report of IEEE Proc. on Field-Programmable Custom Computing Machines (FCCM'06), Napa, California (2006)
【非特許文献１１】O. Maslennikow, V. Lepekha, A. Sergyienko, FPGA Implementation of the Conjugate Gradient Method, LNCS 3911, pp.526-533 (2006)
【非特許文献１２】田原一平, 川口 秀樹, 松岡俊祐，線形計算専用計算機・共役傾斜法マシンの開発に関する研究，平成１９年電気・情報関係学会北海道支部連合大会 (2007), p.147.
【非特許文献１３】D. DuBois, A. DuBois, Sparse Matrix-Vector Multiplicationand Conjugate Gradient Algorithm on Hybrid Computing Platforms, LALP-07-041 (2007)
【発明の開示】
【発明が解決しようとする課題】
【００１１】
行列方程式計算用の専用計算機においては、ソフトウェアプログラムベースの計算機で行う処理と全く同じことをハードウェア化することは困難である。このため、従来は、全体の処理の中でハードウェア化に適し、かつ、ハードウェア化により有効に高速化が図れる処理を適切に選び出してそれらを専用計算機で実行させ、ホスト計算機と連動しながら全体として処理を高速化する方法が採られてきた。
【００１２】
非特許文献４では、密行列をターゲットとしガウスの消去法のための専用計算機を提案している。とりわけ、ガウスの消去法の中で演算量のほとんどを占める前進消去の内積演算の部分をハードウェア化する方式について検討し、計算機のシステム構成を提案している。しかしながら、メモリ容量の制約から、１０００次元以下の行列を対象としており、大規模な計算を行うにはさらなる検討が必要である。
【００１３】
非特許文献５では、同じく密行列をターゲットとしたＬＵ分解法の専用計算機について、ＬＵ分解法の計算ではデータ量（行列の次元）Ｎに対して、演算量がＮ2 、Ｎ1.5 とメモリアクセスがボトルネックになりづらい性質に着目して、ＬＵ分解法に現れる内積演算に特化し並列に高速処理することができる専用計算機の方式検討を行い、計算機のシステム構成を提案している。
【００１４】
また、非特許文献６では、同じく密行列をターゲットとしたＬＵ分解法の専用計算機について、ＬＵ分解法の演算のうち行列のＬＵ分解のための行列成分の四則演算処理のみをハードウェア化し、ピボット選択や前進代入・後退代入などの処理はホスト計算機で処理する方式について検討し、計算機のシステム構成を提案している。
【００１５】
このとき、上記の密行列を対象とした専用計算機による行列計算の高速化の試みにおいては、必然的に大容量のメモリが必要となり、また、より高速計算を実現するためにはハードウェアが大規模となるため、高々、数千次元程度の行列サイズを取り扱うのが限界であり、科学技術計算にしばしば現れるより大規模な計算への適用は困難であるという問題があった。
【００１６】
一方、非特許文献１〜３では、疎行列をターゲットとし疎行列用ＬＵ分解法のための専用計算機を提案している。疎行列性を利用し、行列の次元Ｎに比してはるかに少ない数の非ゼロ行列成分の分だけのローカルメモリを有する演算ユニットを用意し、１行分の加減算・乗除算を並列に行える構成とし、ＬＵ分解法の全ての処理を実行する方式について検討し、計算機のシステム構成を提案している。
【００１７】
また、非特許文献７では、同じく疎行列をターゲットとしたＬＵ分解法の専用計算機について、汎用の算術演算器を２つ用意し、かつ、メモリアクセスを２チャネルとした上でパイプライン処理が可能なハードウェアを小規模ながら実際のＰＣＩインターフェース機能を有するＦＰＧＡ評価キットに実装し、その動作確認まで行っている。
【００１８】
しかしながら、上記のＬＵ分解法をハードウェア化する試みにおいては、ピボット選択も含めた複雑な行列のＬＵ分解、そして、それに続く前進および後退代入の複数の処理が必要であり、ＣＧ法やＢｉＣＧ法などの反復解法に比べて演算量の多い計算を余儀なくされるという問題があった。
【００１９】
非特許文献８では、ラプラス・ポアソン方程式に特化した専用計算機を提案している。この専用計算機では、ラプラス・ポアソン方程式の差分式に現れる大型疎行列の１行ごとの非ゼロ成分の数が高々５であり、かつ、同じ行に現れる行列成分が隣り合うグリッドのアドレスとなることを利用し、この構造に特化した単純反復法のハードウェア化が提案されている。しかし、全体のシステムとしてシンプルなアーキテクチャで実現することができることが判明したものの、やはり、単純反復法は解の収束が悪く、行列解法のＨＰＣ化としては必ずしも適切なものとはなっていなかった。
【００２０】
非特許文献９では、大型疎行列をターゲットとしＣＧ法のための専用計算機を提案している。とりわけ、ＣＧ法の中で行列と探索方向のベクトルとの乗算演算の部分のみをハードウェア化し、他の処理はホスト計算機で実行させ、ホスト計算機と専用計算機とを連動させながらＣＧ法の計算を高速化する方式について検討し、計算機のシステム構成を提案している。
【００２１】
非特許文献１０では、大型疎行列をターゲットとしたＣＧ法の専用計算機について、ＦＰＧＡへの実装を想定した場合の精度とハードウェアサイズ・計算速度のトレードオフについて検討し、パイプライン処理を併用した適切な演算方法を提案している。
【００２２】
非特許文献１１では、大型疎行列をターゲットとしたＣＧ法の専用計算機について、ＣＧ法の演算に現れる除算の効率的な計算方法について検討し、演算部分について、実際にＦＰＧＡに実装しながらその有効性を示している。
【００２３】
非特許文献１２では、大型疎行列をターゲットとしたＣＧ法の専用計算機について、ＣＧ法の１回の反復計算のハードウェア化について検討し、小規模ながら実際にＦＰＧＡに実装しながらその実現可能性を示している。
【００２４】
非特許文献１３では、大型疎行列をターゲットとしたＣＧ法の専用計算機について、商用のリコンフィギュレーション可能な計算機を用いて行列方程式計算のＨＰＣ化について検討し、ＣＧ法に特化したシステムを提案している。
【００２５】
以上のように、従来の行列計算専用計算機は、必ずしも科学技術計算にしばしば現れる大型疎行列には適していないＬＵ分解法をハードウェア化するものや、あるいは大型疎行列に適したＣＧ法のハードウェア化である場合も、ＣＧ法の演算に除算などの複雑な処理が含まれるために、演算方法を部分的にハードウェア化し、ホスト計算機と連動して動作させる構成のものである。これらの行列計算専用計算機は、いずれも、完全な専用計算機という形態ではなく、ホスト計算機とのデータのやり取りを生じるなど動作速度、大規模計算化などの点で必ずしも十分なものではなかった。また、ＣＧ法は基本的に対称行列が対象であるが、科学技術計算では非対称行列もよく現れるため、対称行列に特化した専用計算機では適用範囲が限られてしまうという問題もあった。
【００２６】
そこで、この発明が解決しようとする課題は、電磁界解析、流体解析、構造解析などの各種の数値シミュレーションにおいて最終的に解くべき方程式として現れる大型疎行列方程式の計算を高速に行うことができ、しかもこの大型疎行列方程式の大型疎行列が対称行列である場合および非対称行列である場合の双方に適用することができる行列方程式計算装置および行列方程式計算方法を提供することである。
【課題を解決するための手段】
【００２７】
上記課題を解決するために、第１の発明は、
対象となる大型疎行列方程式を反復解法により解く行列方程式計算装置であって、
上記大型疎行列方程式の大型疎行列の成分のうち、非ゼロの成分のみを１列ごとに格納する所定の個数のメモリを有するメモリ部と、
上記反復解法の演算の少なくとも一部をデータフロー形式で実行する１つまたは複数の演算部とを有し、
上記メモリ部から１行ごとの反復解法の演算に必要なデータを上記演算部に一度にロードするように構成されていることを特徴とするものである。
この行列方程式計算装置は、典型的には、パーソナルコンピュータなどからなるホスト計算機と接続され、このホスト計算機から必要なデータを取得するが、これに限定されるものではなく、必要に応じて、ホスト計算機の機能を持たせてもよい。
【００２８】
第２の発明は、
対象となる大型疎行列方程式を反復解法により解く行列方程式計算方法であって、
上記大型疎行列方程式の大型疎行列の成分のうち、非ゼロの成分のみを１列ごとにメモリ部の所定の個数のメモリに格納するステップと、
上記メモリ部から１行ごとの反復解法の演算に必要なデータを１つまたは複数の演算部に一度にロードするステップと、
上記反復解法の演算の少なくとも一部を上記演算部でデータフロー形式で実行するステップとを有することを特徴とするものである。
【００２９】
第１および第２の発明において、大型疎行列とは、電磁界解析、流体解析、構造解析などの数値シミュレーションにおいて最終的に解くべき行列方程式に現れる数千〜数百万あるいはそれ以上の次元であるがその成分のほとんどはゼロであり、非ゼロ成分は、１行あたり数個〜十数個であるような行列を意味する。大型疎行列は対称行列であっても非対象行列であってもよい。反復解法には、例えば、ＣＧ法、ＢｉＣＧ法、ＩＣＣＧ法、ＢｉＣＧ−Ｓｔａｂ法、ＣＧＳ法、ＢｉＣＧ−Ｓｔａｂ２法などを用いることができる。
【００３０】
第１および第２の発明においては、大型疎行列の成分のうち、非ゼロの成分のみを１列ごとに格納する所定の個数（複数）（必要に応じて選ばれるが、例えば８あるいは１６）のメモリを用いるが、これは次のような理由による。図１は大型疎行列方程式の１例を示す。図１に示すように、行列Ａおよび非同次ベクトルｂからなる大型疎行列方程式
【数１】
から未知数ベクトルｘを求める問題において、これをＣＧ法の計算方法
【数２】
を用いて解く場合を考える。（２−２）式の反復部分の計算には、
【数３】
なる行列とベクトルとの乗算がある。このベクトルｑのｉ番目の成分ｑi の計算には行列Ａのｉ行目の１行分の成分の値が必要であるが、行列の成分の値を１つのメモリに格納しては、必要な値を１度にこれらのメモリから取得することはできない。このため、第１および第２の発明では、メモリ部が、行列の非ゼロ成分のみを１列ごとに格納する複数のメモリを有する構成とし、１行ごとの計算に必要なデータをメモリ部から演算部に並列に一度にロードすることとしている。また、後述の複数の演算部による並列処理の場合は、これらの列ごとにそれぞれ設けられたメモリをさらに並列数分設け、それぞれの並列処理が分担する行の成分の値を格納する。
【００３１】
好適には、行列方程式計算装置は、同じハードウェア回路の複数の演算部を有する。そして、これらの演算部の数だけ、行列の行数Ｎを分割し、それぞれの演算部が連動動作しながら並列処理を行う。こうすることで、行列方程式の計算を演算部の数だけスケーラブルに並列処理により高速化することができる。
好適には、反復解法の演算を可能な限りデータフロー形式で実行する。
【００３２】
さらに、好適には、対称行列を対象としたＣＧ法専用計算機の場合は、偶数かつ複数の同じハードウェア回路の演算部を有する。そして、これらの演算部の２つをペアに連動させることにより、１つの非対称行列を対象としたＢｉＣＧ法専用計算機としても機能させる。また、この演算部のペアの数だけ行列の行数を分割し、それぞれの演算部が連動動作しながら並列処理を行うようにしてもよい。こうすることで、１つの専用計算機を対称行列および非対称行列両方に適用することができ、対称行列モードの場合は非対称行列モードの場合の２倍の速度で動作する計算装置を実現することができる。
【発明の効果】
【００３３】
この発明によれば、行列方程式計算装置にＯＳなどのソフトウェア的な制約なしに、ハードウェア的に許される限り大容量のメモリを搭載可能となるため、大規模な計算を実行することができる。また、演算部におけるＣＧ法などの反復解法のデータフローアーキテクチャ回路構成とメモリ部における行列の１列の非ゼロ成分ごとのメモリ構成とにより、スループットとしても、例えば、ＣＧ法やＢｉＣＧ法では１回の反復計算を３Ｎ＋２クロック（Ｎは行列の次数）で、ＢｉＣＧ−Ｓｔａｂ法やＣＧＳ法では５Ｎクロックで処理することができ、計算方法の限界まで動作時間を最小化した高速な動作を実現することができる。
また、複数の演算部を有する構成とし並列動作させることにより、スケーラブルな並列処理で高速化を図ることができる。
さらに、ＣＧ法の演算部を２つ連動させることでＢｉＣＧ法の演算部を実現することができるため、用途に適した効率的な計算を行うことができる。
さらに、複数の反復解法の演算部を用意し、これらの演算部に同じ行列計算を同時に行わせ、最も早く計算が完了した演算部の解をホスト計算機がアップロードすることにより、自動的に問題に最適な反復解法の方法を選択して行列方程式の計算を実行することができる。
【発明を実施するための最良の形態】
【００３４】
以下、この発明の実施の形態について図面を参照しながら説明する。
まず、この発明の第１の実施の形態による行列方程式計算装置について説明する。
この行列方程式計算装置はＣＧ法専用計算機により構成される。この行列方程式計算装置は、（１）式の行列方程式を（２−１）式、（２−２）式および（３）式に示すＣＧ法により計算するものである。
【００３５】
図２はこの行列方程式計算装置１の全体構成、図３はこの行列方程式計算装置１の各部の構成の詳細を示す。図２および図３に示すように、この行列方程式計算装置１は、大型疎行列Ａの全成分の値を格納するメモリモジュール１１と、（２−２）式のＣＧ法の１回分の反復計算をデータフロー形式で３Ｎクロックで実行する演算モジュール１２と、メモリモジュール１１と演算モジュール１２との間で、（３）式の大型疎行列Ａとベクトルｐk との乗算をクロック遅延のない組合せ回路で実行し、この結果の値Ａｐk を演算モジュール１２に提供する行列−ベクトル乗算回路１３と、演算モジュール１２で得られたｋ＋１番目の反復解が収束しており最終的な解となっているかどうかを判定する収束判定部１４とを有する。メモリモジュール１１、行列−ベクトル乗算回路１３、演算モジュール１２および収束判定部１４における一連のデータの流れの動作はマスターコントローラ１５により制御される。
【００３６】
メモリモジュール１１は、大型疎行列Ａを非ゼロ成分のみに圧縮したＮ×ｍのサイズの行列Ａ’の成分を列ごとに格納するｍ個のメモリＭ１−１〜Ｍ１−ｍを有する。この大型疎行列Ａの圧縮は具体的には次のように行う。例えば、各行ごとの非ゼロ成分の数のうち全Ｎ行の中で最大のものをｍとしてＮ×ｍのサイズの行列を用意し、これに各行ごとに非ゼロ成分を左詰めで格納する。同時に、これらのメモリＭ１−１〜Ｍ１−ｍと全く同じサイズのメモリをもう１つずつ用意する（メモリＭ２−１〜Ｍ２−ｍ）。これらのメモリＭ２−１〜Ｍ２−ｍには、対応するメモリＭ１−１〜Ｍ１−ｍの非ゼロ成分がもともとの大型疎行列Ａのどの列に位置していたかを示すインデックスを格納する。これにより、通常、大型疎行列ではＮが数千〜数千万に対しｍは高々十数個であるため、大型疎行列Ａの情報を失うことなく大幅にメモリ容量を節約することができる。メモリモジュール１１をこのように構成することで、メモリ容量を大幅に節約することができ、かつ、一度に１行分の行列の非ゼロ成分全てに並列アクセスすることができるためスループットとしての計算性能に大きく影響する行列の成分が格納されたメモリＭ１−１〜Ｍ１−ｍへのアクセスは、最小限の１回にまで減らすことができる。
【００３７】
メモリモジュール１１はさらに、ＣＧ法の反復計算実行時に一時的に値を格納しておくために、同じベクトルｐk の値をｍ個別々に格納するためのメモリＭ３−０〜Ｍ３−ｍ、ベクトルｐk の更新値ｐk+1 を格納するためのメモリＭ３−０’、ベクトルｘk を格納するためのメモリＭ４、ベクトルｘk の更新値ｘk+1 を格納するためのメモリＭ４’、ベクトルｒk を格納するためのメモリＭ５、ベクトルｒk の更新値ｒk+1 を格納するためのメモリＭ５’、ベクトルｑk を格納するためのメモリＭ６、ベクトルｒk のノルムの２乗‖ｒk ‖2 ＝（ｒk ，ｒk ）を格納するためのメモリＭ７、ノルムの２乗‖ｒk ‖2 の更新値‖ｒk+1 ‖2 ＝（ｒk+1 ，ｒk+1 ）を格納するためのメモリＭ７’、行列の次数Ｎ、行ごとの非ゼロ成分の最大数ｍ、非同次項ベクトルｂのノルムの２乗‖ｂ‖2 ＝（ｂ，ｂ）、相対残差の収束判定基準値ε、最大反復回数Ｋｍ、各反復ステップの残差Ｒｓなどの反復の収束条件を格納するためのメモリＭ８を有する。
【００３８】
上述のメモリＭ１−１〜Ｍ１−ｍ、Ｍ２−１〜Ｍ２−ｍ、Ｍ３−０〜Ｍ３−ｍ、Ｍ３−０’、Ｍ４、Ｍ４’、Ｍ５、Ｍ５’、Ｍ６、Ｍ７、Ｍ７’、Ｍ８としては、それぞれ独立した個別のメモリを用いてもよいし、それらのうちの１つまたは２つ以上のメモリとして、独立した１つのメモリのメモリ領域を所定個数に分割したものを用い、これらのメモリ領域をメモリＭ１−１〜Ｍ１−ｍ、Ｍ２−１〜Ｍ２−ｍ、Ｍ３−０〜Ｍ３−ｍ、Ｍ３−０’、Ｍ４、Ｍ４’、Ｍ５、Ｍ５’、Ｍ６、Ｍ７、Ｍ７’、Ｍ８のうちのいずれかに用いてもよい。これらのメモリＭ１−１〜Ｍ１−ｍ、Ｍ２−１〜Ｍ２−ｍ、Ｍ３−０〜Ｍ３−ｍ、Ｍ３−０’、Ｍ４、Ｍ４’、Ｍ５、Ｍ５’、Ｍ６、Ｍ７、Ｍ７’、Ｍ８としては、例えば、ＳＲＡＭまたはＤＲＡＭを用いることができる。
【００３９】
行列−ベクトル乗算回路１３は、上述のように列ごとにメモリＭ１−１〜Ｍ１−ｍに分割格納された大型疎行列Ａの成分を、メモリＭ２−１〜Ｍ２−ｍを参照しつつ１行ごとに１行の全ての成分に並列に同時アクセスし、また、ベクトルｐk を格納するメモリＭ３−１〜Ｍ３−ｍからそれぞれ（３）式の計算に必要な値のアドレスに同時アクセスし、データフロー形式で（３）式の計算を実行し、ベクトルｑk を１クロックで１成分ずつ計算し、その結果を演算モジュール１２に引き渡す。
【００４０】
演算モジュール１２は、行列−ベクトル乗算回路１３からのベクトルｑk 、メモリＭ３−０のベクトルｐk 、メモリＭ４のベクトルｘk 、メモリＭ５のベクトルｒk 、メモリＭ６のベクトルｑk およびメモリＭ７の（ｒk ，ｒk ）を用いて、（２−２）式（ｉ）中のベクトルｐk とベクトルｑk （＝Ａｐk ）との内積（ｐk ，ｑk ）の計算を実行するための演算回路３１、（２−２）式（ｉ）中の除算（ｒk ，ｒk ）／（ｐk ，ｑk ）の計算を実行するための演算回路３２、（２−２）式（ｉｉ）のベクトルｘk+1 の計算を実行するための演算回路３３、（２−２）式（ｉｉｉ）のベクトルｒk+1 の計算を実行するための演算回路３４、（２−２）式（ｖ）中のベクトルｒk+1 の２乗（ｒk+1 ，ｒk+1 ）の計算を実行するための演算回路３５、（２−２）式（ｉｖ）中の除算（ｒk+1 ，ｒk+1 ）／（ｒk ，ｒk ）の計算を実行するための演算回路３６、（２−２）式（ｖ）のベクトルｐk+1 の計算を実行するための演算回路３７を有する。
【００４１】
（２−２）式からわかるように、演算回路３１、３３、３４、３５、３７の計算はそれぞれＮクロックで、演算回路３２、３６の計算は１クロックで実行することができるが、演算モジュール１２は全体として可能な限りこれらの計算を並列に実行することができるように構成されている。具体的には、演算回路３１、３２のＮ＋１クロックの計算が完了した後は、演算回路３３、３４、３５の計算は、演算回路３４、３５のデータフロー性を利用するとまとめてＮクロックで実行することができるため、残りの演算回路３６、３７のＮ＋１クロックの計算を合わせて全体で３Ｎ＋２クロックで実行することができるように構成されている。
【００４２】
収束判別部１４は、ＣＧ法の演算を実行する行列−ベクトル乗算回路１３から演算モジュール１２での１回分の反復計算で得られたｋ＋１回目の残差ベクトルｒk+1 から相対誤差を計算し、あらかじめ設定されたメモリＭ８内の相対残差の収束判定基準値ε内に収まり反復計算が収束しているかどうかを判別する。
【００４３】
マスターコントローラ１５は、計算が収束していればｘk+1 を最終的な近似解として計算を終了し、ホスト計算機（図示せず）に正常に収束値を得た旨を通知する。計算が収束していなければ、あらかじめ設定されたメモリＭ８内の最大反復回数Ｋｍを超えた場合は、ホスト計算機に反復計算が収束しなかった旨を通知し、最大反復回数Ｋｍを超えていない場合は、反復計算で得られたベクトルｐk+1 、ｘk+1 、ｒk+1 および（ｒk+1 ，ｒk+1 ）を次の反復計算に送り、再度、行列−ベクトル乗算回路１３から演算モジュール１２での反復計算の実行を指示する。
【００４４】
図４は行列方程式計算装置１の使用形態を示す。行列方程式計算装置１はＰＣＩインターフェース４１によりホスト計算機４２と接続されている。ＣＧ法の演算を実行する上で必要な値は、圧縮格納された大型疎行列Ａ、反復計算の初期値となるベクトルｐ0 （＝ｒ0 ）、ｘ0 、非同次項の（ｂ，ｂ）、収束判定基準ε、最大反復回数Ｋｍである。これらの値はホスト計算機４２に格納され、計算開始指示の前に行列方程式計算装置１にダウンロードされてメモリモジュール１１のメモリＭ３、Ｍ４、Ｍ５、Ｍ８に格納されるようになっている。ホスト計算機４２は、これらの値を行列方程式計算装置１にダウンロードした後、行列方程式計算装置１に対し計算開始を指示する。行列方程式計算装置１は、（２−２）式の反復計算および収束判定を繰返し、最大反復回数Ｋｍ以内で解が収束判定基準ε以下になれば、メモリＭ４’のベクトルｘk+1 を近似解としてホスト計算機４２に引渡し、最大反復回数Ｋｍまで計算しても収束しなければ、計算は未収束として、これをホスト計算機４２に通知する。
【００４５】
図５Ａ〜Ｅは行列方程式計算装置１の実装例を示す。図５ＡおよびＢに示すように、この例では、１つの小プリント基板５１上に、メモリモジュール１１の２列分の圧縮行列の成分を格納するメモリ５２と、この２列の行列と対応する探索ベクトルｐk の成分との乗算を行う行列−ベクトル乗算回路１３とが搭載される。図５ＣおよびＤに示すように、この小プリント基板５１が合計１６個、サブプリント基板上５３に接続される。図５Ｅに示すように、メインプリント基板上５４に、演算モジュール１２、収束判定部１４およびマスターコントローラ１５が搭載され、これにサブプリント基板５３が接続される。こうして、ＣＧ法専用計算機としての行列方程式計算装置１が構成される。これらの搭載部品のほとんどは例えばＦＰＧＡまたはＡＳＩＣにより作製することができる。このように構成することにより、行列方程式計算装置１をある決められた列数ｍ用に構成した場合でも、あらかじめ十分なソケットを用意しておけば、後に小プリント基板５１を追加接続することにより、行列サイズの変更に対してフレキシビリティーを持たすことができる。
【００４６】
この第１の実施の形態によれば、次のような利点を得ることができる。すなわち、行列方程式計算装置１をＣＧ法専用計算機により構成し、メモリモジュール１１と演算モジュール１２とを分離した構成としていることにより、数値シミュレーションに用いる大型疎行列に応じてメモリモジュール１１を適切な構成のものに交換することができ、比較的複雑な構成を有する演算モジュール１３は変更しないで、大型疎行列の変更に柔軟に対応することができ、フレキシビリティーが高い。また、メモリモジュール１１は、行列の１列ごとにｍ個のメモリＭ１−１〜Ｍ１−ｍを設け、これらのメモリＭ１−１〜Ｍ１−ｍから（３）式の１行分の計算に必要なデータを１クロックで行列−ベクトル乗算回路１３にロードし、（３）式の１行分の計算を並列的に実行するようにしていることにより、メモリモジュール１１へのアクセス回数を最小限にすることができ、それによって高速動作化を図ることができ、行列方程式の計算を高速かつ短時間で行うことができる。また、この行列方程式計算装置１は実装や計算の大規模化が容易である。以上により、極めて実用性が高い高性能の行列方程式計算装置１を実現することができる。
この行列方程式計算装置１は、電磁界解析、流体解析、構造解析などの各種の数値シミュレーションで最終的に解くべき方程式として現れる大型疎行列方程式を解くのに適用して好適なものである。
【００４７】
次に、この発明の第２の実施の形態による行列方程式計算装置について説明する。
この第２の実施の形態においては、メモリモジュール１１へのアクセス回数の最小限化に加え、さらに、演算モジュール１２における３Ｎ＋２クロック分の計算をより高速化すべく、演算モジュール１２を複数用意し、３Ｎ＋２クロック分の計算をこれらの複数の演算モジュール１２で並列に実行する方式が用いられる。演算モジュール１２における３Ｎ＋２クロックの計算はもともと全て並列処理してもよく、この部分の計算を複数の演算モジュール１２で行うことにより、１つの演算モジュール１２を用いた場合に比べ、演算モジュール１２の台数分だけ、高速化を図ることができる。上記以外のことは第１の実施の形態と同様である。
この行列方程式計算装置１の実装は、図６に示すように、図５Ｅに示すものと同様なメインプリント基板５４上にサブプリント基板５３を複数搭載することにより簡単に行うことができる。
【００４８】
この第２の実施の形態によれば、第１の実施の形態と同様な利点に加えて、動作速度のより一層の向上を図ることができ、行列方程式の計算をより高速に行うことができるという利点を得ることができる。
【００４９】
次に、この発明の第３の実施の形態による行列方程式計算装置について説明する。
この第３の実施の形態においては、複数の演算モジュール３の並列処理による高速化に加え、さらに、ＣＧ法は大型疎行列が対称行列である場合のみに適用が限定されているのに対し、非対称行列である場合にも適用することができる機能を装備すべく、２つの演算モジュール１２を連動動作させ、非対称行列にも適用することができるＢｉＣＧ法マシンとして動作させる方式が用いられる。ＢｉＣＧ法はもともとＣＧ法に転置行列に対するＣＧ法の処理を追加したものであり、もう１つの演算モジュール１２に転置行列を処理させることにより、２つの演算モジュール１２をペアで連動動作させることによりＢｉＣＧ法の動作を実現することができる。上記以外のことは第１の実施の形態と同様である。
この行列方程式計算装置１の実装は、図６に示すように、図５Ｅに示すものと同様なメインプリント基板５４上にサブプリント基板５３を複数搭載することにより簡単に行うことができる。
【００５０】
次に、ハードウェア記述言語であるＶＨＤＬ（VHSIC Hardware Description Language)により行列方程式計算装置１による大型疎行列計算の論理シミュレーションを行った結果について説明する。
図７はシミュレーションの例として使用したリニアモータの界磁部のモデルを、図８はその数値モデルを示す。このモデルは、３つの励磁コイル６１〜６３で励磁された静磁場が、鉄（μr ＝５０００）からなる磁場遮蔽材料６４によって外部にもれないよう遮蔽される現象に関するもので、系は軸対称性を持っており、図８の数値シミュレーションの解析領域は軸対称性を考慮して２次元の５０×１００グリッドサイズのものである。この磁場を記述する支配方程式は、与えられた電流密度ベクトルＪ、透磁率μに対し、ベクトルポテンシャルＡを未知数とし、
【数４】
となる。この方程式を半径方向、軸方向ともにΔｌのサイズの一様なグリッドで分割した上で軸対称性を考慮して差分法に基づき差分式で表すと、半径方向にi 番目、軸方向にj
番目のグリッドでは、
【数５】
となる。ただし、μi,j 、Ａi,j 、Ｊi,j はそれぞれ、半径方向にi 番目、軸方向にj 番目のグリッドでの透磁率μ、ベクトルポテンシャルＡのθ成分、電流密度ベクトルＪのθ成分である。このとき、未知数Ａi,j の数は、全てのグリッド数の５０×１００＝５０００であるので、差分法により数値解析する際に現れる大型疎行列は５０００×５０００のサイズとなる。この大型疎行列計算を、非対称行列用のＢｉＣＧ法に基づいて構成した行列方程式計算装置１により実行する場合のＶＨＤＬ論理回路シミュレーションにより行った。その解から得られた静磁場分布の結果を図９に示す。実際の計算では、（５）式からわかるように、１行に現れる未知数Ａi,j の係数は高々５個であるので、この疎行列性を利用し、行列は５０００×５のサイズに圧縮したものを用いている。図９に示す結果は、同じ数値モデルに基づいてＣ言語ソフトウェアシミュレーションを行った結果と一致しており、この行列方程式計算装置１の妥当性が保証される。このＶＨＤＬ論理回路シミュレーションでの動作を参考に、この行列方程式計算装置１がアクセス速度２６６ＭＨｚのメモリと同程度の速度で動作した場合の性能を見積ると、２．６ＧＨｚのＣＰＵを有するパソコン上でＣ言語を用いて計算した時の約５０〜６０倍の性能に匹敵する。
【００５１】
ここでは、簡単のため、軸対称２次元となるモデルを用いたが、基本的には３次元の場合も全く同様である。すなわち、例えば同じ図７の３次元モデルでは、グリッドサイズが１００×１００×１００であるので、最終的に現れる大型疎行列は、１００００００×１００００００で、（４）式に対応する３次元の差分式では、係数は１６個となるため、最終的には、１００００００×１６のサイズに圧縮された行列を解くことに帰着する。このため、この行列方程式計算装置１では、軸対称２次元の場合が、Ｎ＝５０００、ｍ＝５であったものが、３次元では、単に、Ｎ＝１００００００、ｍ＝１６となるだけで、動作自体は全く同様である。
この第３の実施の形態によれば、第１の実施の形態と同様な利点に加えて、大型疎行列Ａが非対称行列である場合にも行列方程式計算装置１を適用することができるという利点を得ることができる。
【００５２】
以上、この発明の実施の形態について具体的に説明したが、この発明は、上述の実施の形態に限定されるものではなく、この発明の技術的思想に基づく各種の変形が可能である。
例えば、上述の実施の形態において挙げた数値、回路構成、配置などはあくまでも例に過ぎず、必要に応じて、これらと異なる数値、回路構成、配置などを用いてもよい。また、行列方程式の解法には、ＣＧ法やＢｉＣＧ法ではなく、ＩＣＣＧ法、ＢｉＣＧ−Ｓｔａｂ法、ＣＧＳ法、ＢｉＣＧ−Ｓｔａｂ２法などの他の反復解法を採用し、演算モジュール１２もこれらの行列方程式の解法用のもので構成してもよい。
【図面の簡単な説明】
【００５３】
【図１】大型疎行列方程式の１例を示す略線図である。
【図２】この発明の第１の実施の形態による行列方程式計算装置の全体構成を示す略線図である。
【図３】この発明の第１の実施の形態による行列方程式計算装置の各部の構成の詳細を示す略線図である。
【図４】この発明の第１の実施の形態による行列方程式計算装置をホスト計算機と接続した状態を示す略線図である。
【図５】この発明の第１の実施の形態による行列方程式計算装置の実装例を示す略線図である。
【図６】この発明の第２の実施の形態による行列方程式計算装置の実装例を示す略線図である。
【図７】この発明の第３の実施の形態による行列方程式計算装置を用いて行った数値シミュレーションに用いた数値モデルを示す略線図である。
【図８】この発明の第３の実施の形態による行列方程式計算装置を用いて行った数値シミュレーションに用いた数値モデルを示す略線図である。
【図９】この発明の第３の実施形態による行列方程式計算装置を用いて行った数値シミュレーションの結果を示す略線図である。
【図１０】電磁界解析を例に、多くの科学技術計算は最終的に大型疎行列方程式の計算に帰着することを説明するための略線図である。
【符号の説明】
【００５４】
１…行列方程式計算装置、１１…メモリモジュール、１２…演算モジュール、１３…行列−ベクトル乗算回路、１４…収束判定部、１５…マスターコントローラ、４２…ホスト計算機、５１…プリント基板、５３…サブプリント基板、５４…メインプリント基板、Ｍ１−１〜Ｍ１−ｍ、Ｍ２−１〜Ｍ２−ｍ、Ｍ３−０〜Ｍ３−ｍ、Ｍ３−０’、Ｍ４、Ｍ４’、Ｍ５、Ｍ５’、Ｍ６、Ｍ７、Ｍ７’、Ｍ８…メモリ
【特許請求の範囲】
【請求項１】
対象となる大型疎行列方程式を反復解法により解く行列方程式計算装置であって、
上記大型疎行列方程式の大型疎行列の成分のうち、非ゼロの成分のみを１列ごとに格納する所定の個数のメモリを有するメモリ部と、
上記反復解法の演算の少なくとも一部をデータフロー形式で実行する１つまたは複数の演算部とを有し、
上記メモリ部から１行ごとの反復解法の演算に必要なデータを上記演算部に一度にロードするように構成されていることを特徴とする行列方程式計算装置。
【請求項２】
上記演算部を複数有することを特徴とする請求項１記載の行列方程式計算装置。
【請求項３】
上記複数の演算部を並列動作させることを特徴とする請求項２記載の行列方程式計算装置。
【請求項４】
上記反復解法は共役傾斜法または双共役傾斜法であることを特徴とする請求項３記載の行列方程式計算装置。
【請求項５】
対象となる大型疎行列方程式を反復解法により解く行列方程式計算方法であって、
上記大型疎行列方程式の大型疎行列の成分のうち、非ゼロの成分のみを１列ごとにメモリ部の所定の個数のメモリに格納するステップと、
上記メモリ部から１行ごとの反復解法の演算に必要なデータを１つまたは複数の演算部に一度にロードするステップと、
上記反復解法の演算の少なくとも一部を上記演算部でデータフロー形式で実行するステップとを有することを特徴とする行列方程式計算方法。
【請求項６】
上記演算部を複数用いることを特徴とする請求項５記載の行列方程式計算方法。
【請求項７】
上記複数の演算部を並列動作させることを特徴とする請求項６記載の行列方程式計算方法。
【請求項８】
上記反復解法は共役傾斜法または双共役傾斜法であることを特徴とする請求項７記載の行列方程式計算方法。
【請求項１】
対象となる大型疎行列方程式を反復解法により解く行列方程式計算装置であって、
上記大型疎行列方程式の大型疎行列の成分のうち、非ゼロの成分のみを１列ごとに格納する所定の個数のメモリを有するメモリ部と、
上記反復解法の演算の少なくとも一部をデータフロー形式で実行する１つまたは複数の演算部とを有し、
上記メモリ部から１行ごとの反復解法の演算に必要なデータを上記演算部に一度にロードするように構成されていることを特徴とする行列方程式計算装置。
【請求項２】
上記演算部を複数有することを特徴とする請求項１記載の行列方程式計算装置。
【請求項３】
上記複数の演算部を並列動作させることを特徴とする請求項２記載の行列方程式計算装置。
【請求項４】
上記反復解法は共役傾斜法または双共役傾斜法であることを特徴とする請求項３記載の行列方程式計算装置。
【請求項５】
対象となる大型疎行列方程式を反復解法により解く行列方程式計算方法であって、
上記大型疎行列方程式の大型疎行列の成分のうち、非ゼロの成分のみを１列ごとにメモリ部の所定の個数のメモリに格納するステップと、
上記メモリ部から１行ごとの反復解法の演算に必要なデータを１つまたは複数の演算部に一度にロードするステップと、
上記反復解法の演算の少なくとも一部を上記演算部でデータフロー形式で実行するステップとを有することを特徴とする行列方程式計算方法。
【請求項６】
上記演算部を複数用いることを特徴とする請求項５記載の行列方程式計算方法。
【請求項７】
上記複数の演算部を並列動作させることを特徴とする請求項６記載の行列方程式計算方法。
【請求項８】
上記反復解法は共役傾斜法または双共役傾斜法であることを特徴とする請求項７記載の行列方程式計算方法。
【図１】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【図１０】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【図１０】
【公開番号】特開２０１０−１２２８５０（Ｐ２０１０−１２２８５０Ａ）
【公開日】平成２２年６月３日（２０１０．６．３）
【国際特許分類】
【出願番号】特願２００８−２９５１３５（Ｐ２００８−２９５１３５）
【出願日】平成２０年１１月１９日（２００８．１１．１９）
【出願人】（５０４１９３８３７）国立大学法人室蘭工業大学 (70)
【Ｆターム（参考）】
【公開日】平成２２年６月３日（２０１０．６．３）
【国際特許分類】
【出願日】平成２０年１１月１９日（２００８．１１．１９）
【出願人】（５０４１９３８３７）国立大学法人室蘭工業大学 (70)
【Ｆターム（参考）】
[ Back to top ]