JP3809062B2 - マルチレベル不完全ブロック分解による前処理を行う処理装置 - Google Patents
マルチレベル不完全ブロック分解による前処理を行う処理装置 Download PDFInfo
- Publication number
- JP3809062B2 JP3809062B2 JP2000354179A JP2000354179A JP3809062B2 JP 3809062 B2 JP3809062 B2 JP 3809062B2 JP 2000354179 A JP2000354179 A JP 2000354179A JP 2000354179 A JP2000354179 A JP 2000354179A JP 3809062 B2 JP3809062 B2 JP 3809062B2
- Authority
- JP
- Japan
- Prior art keywords
- matrix
- coefficient
- variable
- storage means
- elements
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Expired - Fee Related
Links
- 238000000354 decomposition reaction Methods 0.000 title claims description 54
- 238000007781 pre-processing Methods 0.000 title claims description 28
- 239000011159 matrix material Substances 0.000 claims description 221
- 238000000034 method Methods 0.000 claims description 72
- 239000013598 vector Substances 0.000 claims description 40
- 230000008569 process Effects 0.000 claims description 20
- 230000008901 benefit Effects 0.000 claims description 6
- 238000004364 calculation method Methods 0.000 description 22
- 238000006243 chemical reaction Methods 0.000 description 16
- 238000010586 diagram Methods 0.000 description 8
- 230000010365 information processing Effects 0.000 description 7
- 238000004422 calculation algorithm Methods 0.000 description 5
- 238000004458 analytical method Methods 0.000 description 3
- KHOITXIGCFIULA-UHFFFAOYSA-N Alophen Chemical compound C1=CC(OC(=O)C)=CC=C1C(C=1N=CC=CC=1)C1=CC=C(OC(C)=O)C=C1 KHOITXIGCFIULA-UHFFFAOYSA-N 0.000 description 2
- 238000004891 communication Methods 0.000 description 2
- 230000000694 effects Effects 0.000 description 2
- 239000000284 extract Substances 0.000 description 2
- 230000006870 function Effects 0.000 description 2
- 229940050561 matrix product Drugs 0.000 description 2
- 230000003287 optical effect Effects 0.000 description 2
- 238000003491 array Methods 0.000 description 1
- 230000004888 barrier function Effects 0.000 description 1
- 230000005540 biological transmission Effects 0.000 description 1
- 230000000295 complement effect Effects 0.000 description 1
- 230000008030 elimination Effects 0.000 description 1
- 238000003379 elimination reaction Methods 0.000 description 1
- 239000012530 fluid Substances 0.000 description 1
- 230000001737 promoting effect Effects 0.000 description 1
- 230000008707 rearrangement Effects 0.000 description 1
- 230000009467 reduction Effects 0.000 description 1
- 239000004065 semiconductor Substances 0.000 description 1
- 238000012916 structural analysis Methods 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
- G06F17/12—Simultaneous equations, e.g. systems of linear equations
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Mathematical Physics (AREA)
- Computational Mathematics (AREA)
- Mathematical Analysis (AREA)
- Mathematical Optimization (AREA)
- Pure & Applied Mathematics (AREA)
- Data Mining & Analysis (AREA)
- Theoretical Computer Science (AREA)
- Operations Research (AREA)
- Databases & Information Systems (AREA)
- Software Systems (AREA)
- General Engineering & Computer Science (AREA)
- Algebra (AREA)
- Complex Calculations (AREA)
Description
【０００１】
【発明の属する技術分野】
本発明は、連立一次方程式の解を求める反復法に係り、さらに詳しくは、反復法における前処理（preconditioning ）を行う処理装置に関する。
【０００２】
【従来の技術】
物理現象の解析においては、しばしば、大規模な行列を係数行列とする次のような連立一次方程式を解く必要が生じる。
Ａｘ＝ｂ （１）
ここで、Ａはｎ×ｎの係数行列、ｘはｎ次元の変数ベクトル、ｂはｎ次元の定数ベクトルである。ｎの値が１０6 程度になることも珍しくない。
【０００３】
このような大規模な連立一次方程式は、気象予測、原子炉設計、半導体の回路解析、航空工学における流体解析、構造物の構造解析等の多くの科学技術計算に用いられる。したがって、大規模な連立一次方程式を効率よく高速に解くことは、科学技術計算の重要な問題の１つである。
【０００４】
コンピュータを用いて（１）式を解く１つの方法として、ＡをＬＵ分解するガウスの消去法に基づいた直接法がある。しかし、Ａが大きなスパース行列の場合、非零要素が各行に数個しかないこともあり、計算コストや記憶領域の面で無駄が多い。そこで、このような場合には、単純な行列ベクトル積を繰り返して近似解を求める反復法が多く用いられている。
【０００５】
スパース行列の反復法の収束を促進する方法として、前処理が用いられる。この処理では、適当な前処理行列（preconditioner）Ｍ-1を用いて、（１）式が次のように変更される。
Ｍ-1Ａｘ＝Ｍ-1ｂ （２）
そして、（２）式に対して反復法を適用することで、（１）式の近似解が求められる。このとき、ＭはＡの不完全ブロック分解（block incomplete factorization，ＢＩＦ）に対応し、Ｍ-1はＭの逆行列を表す。このように、Ａの不完全ブロック分解を利用して前処理を行う１つの方法として、代数的マルチレベル反復法（algebraic multilevel iteration method ）が知られている。
【０００６】
【発明が解決しようとする課題】
しかしながら、上述した従来の代数的マルチレベル反復法による前処理には、次のような問題がある。
【０００７】
一般に、代数的マルチレベル反復法による前処理は、係数行列ＡがＭ行列（M-matrix）であるような問題に対して適用でき、反復法の収束を速める効果があることが知られている。ここで、ＡがＭ行列であるとは、Ａ＝（ａij）が正則行列であり、さらに次の２つの条件を満たすことを意味する。
（１）ａij≦０（ｉ≠ｊ）（Ａの非対角要素が０もしくは負）
（２）Ａ-1≧０（Ａの逆行列Ａ-1の全要素が非負）
しかし、ＡがＭ行列でない場合は、代数的マルチレベル反復法による前処理はあまり効果がなく、計算が収束しないという問題がある。
【０００８】
本発明の課題は、係数行列がＭ行列でないような連立一次方程式の反復法において、計算の収束を速めるような前処理を行う処理装置を提供することである。
【０００９】
図１は、本発明の処理装置の原理図である。図１の処理装置は、係数格納手段１、決定手段２、集合格納手段３、逆行列手段４、係数行列手段５、分解手段６、および計算手段７を備え、係数行列がＭ行列でない連立一次方程式を解く反復法のための前処理を、マルチレベル不完全ブロック分解により行なう。
【００１０】
係数格納手段１は、不完全ブロック分解のあるレベルにおける係数行列を格納する。決定手段２は、係数格納手段１に格納された係数行列の要素のうち、取り除かれる変数の変数番号を行番号および列番号とする要素からなるブロック行列が対角優位（diagonal dominant ）になるように、それらの取り除かれる変数の変数番号の集合を決定する。このとき、決定手段２は、係数格納手段１に格納された係数行列の要素のうち、対角優位な各要素の添え字に乱数を割り当て、２つの添え字に割り当てられた乱数を比較して大きい方に対応する添え字を取り除く処理を繰り返すことで、残された添え字からなる最大独立集合を求め、その最大独立集合を変数番号の集合とする。そして、集合格納手段３は、決定された変数番号の集合に関する情報を格納する。
【００１１】
また、逆行列手段４は、集合格納手段３に格納された情報に基づいて、上記ブロック行列の近似逆行列を求める。係数行列手段５は、係数格納手段１に格納された係数行列の要素のうち、残される変数の変数番号を行番号および列番号とする要素からなるブロック行列から、次のレベルの係数行列を求める。そして、分解手段６は、次の不完全ブロック分解を行うために、次のレベルの係数行列を係数格納手段１に格納する。
【００１２】
また、逆行列手段４は、集合格納手段３に格納された情報に基づいて、上記ブロック行列の近似逆行列を求める。係数行列手段５は、係数格納手段１に格納された係数行列の要素のうち、残される変数の変数番号を行番号および列番号とする要素からなるブロック行列から、次のレベルの係数行列を求める。そして、分解手段６は、次の不完全ブロック分解を行うために、次のレベルの係数行列を係数格納手段１に格納する。
【００１３】
次に、逆行列手段４は、その集合の情報に基づいて、取り除かれる変数の係数からなるブロック行列の逆行列を近似的に求め、係数行列手段５は、上記係数行列の情報を参照して、処理対象として残される変数の係数からなるブロック行列から、次のレベルの係数行列を計算する。そして、分解手段６が、得られた係数行列を分解される係数行列として係数格納手段１に格納することで、次のレベルの不完全ブロック分解が再帰的に行われる。
【００１４】
計算手段７は、再帰的な不完全ブロック分解により得られた各レベルの近似逆行列を用いて、反復法で必要となる行列ベクトル積を再帰的に計算することで、反復法のための前処理を行う。
【００１５】
対角優位な行列の場合、後述するように、与えられた近似精度で反復計算により、その逆行列を求めることができる。したがって、取り除かれる変数の係数からなるブロック行列が対角優位になるように、それらの変数の変数番号の集合を決定することで、任意のスパース行列を適切な近似精度で分解することができる。このような不完全ブロック分解によれば、Ｍ行列でない係数行列に対しても、反復法の収束を速める前処理が可能になる。
【００１６】
例えば、図１の係数格納手段１および集合格納手段３は、後述する図６の共有メモリ１４または図１０のメモリ３２に対応し、図１の決定手段２、逆行列手段４、係数行列手段５、分解手段６、および計算手段７は、図６のプロセッサ１１または図１０のＣＰＵ（中央処理装置）３１に対応する。
【００１７】
また、取り除かれる変数の変数番号の集合は、後述する集合Ｆに対応し、それらの変数番号を行番号および列番号とする要素からなるブロック行列は、後述する（３）式の行列ＡFFに対応する。また、近似逆行列は、後述する（４）式のＭFF -1または図４のＦFFに対応し、残される変数の変数番号を行番号および列番号とする要素からなるブロック行列は、（３）式の行列ＡCCに対応し、次のレベルの係数行列は、図４の行列ＡC または（４）式のＭC に対応する。
【００１８】
【発明の実施の形態】
以下、図面を参照しながら、本発明の実施の形態を詳細に説明する。
一般に、前処理を施した反復法においては、（２）式のＭ-1とｎ次元ベクトルｑの行列ベクトル積Ｍ-1ｑを繰り返し計算する必要が生じる。本実施形態の処理装置は、代数的マルチレベル反復法による前処理において、連立一次方程式の係数行列のマルチレベル不完全ブロック分解を行い、各レベルで得られた結果を利用してＭ-1ｑを計算する。
【００１９】
まず、連立方程式の変数番号の集合｛１，．．．，ｎ｝から適当に選択された変数番号の集合をＣとし、その補集合をＦとすると、ｎ×ｎの行列Ａは以下の形に変換できる。
【００２０】
【数１】
【００２１】
ここで、集合Ｃは、残される変数を表す変数番号の集合に対応し、集合Ｆは、取り除かれる変数を表す変数番号の集合に対応する。また、ブロック行列ＡCCの要素の行番号および列番号は、集合Ｃの要素により表され、ＡCCの行および列の数は、集合Ｃの要素の数に一致する。同様に、ブロック行列ＡFFの要素の行番号および列番号は、集合Ｆの要素により表され、ＡFFの行および列の数は、集合Ｆの要素の数に一致する。
【００２２】
また、ＡFFの近似逆行列をＭFF -1とし、ＡC ：＝ＡCC−ＡCFＭFF -1ＡFCとすると、ＡC は、例えば、ＢC ＝ＡC −ＲC により近似することができる。ただし、ＲC は、近似により落とされる非ゼロ要素を含む行列である。ＡC の近似をＭC と書くことにすると、Ａの不完全ブロック分解Ｍは以下のように書ける。
【００２３】
【数２】
【００２４】
そして、ＭC の不完全ブロック分解を、粗いレベルからさらに粗いレベルへと再帰的に繰り返すことにより、各レベルにおいて、（４）式のような不完全ブロック分解を得ることができる。また、ｐ：＝Ｍ-1ｑ＝ＢＩＦ（０，ｑ）は、（４）式のブロック行列を用いて再帰的に計算することができる。最初のレベルを０とし、最も粗いレベルをＬとすると、この行列ベクトル積の計算アルゴリズムは、以下のようになる。
このとき、Ｉ＝ＭC -1ＡC （Ｉは単位行列）の近似の度合いは、｜Ｉ−ＭFF -1ＡFF｜、｜ＭC -1ＲC ｜、および｜Ｉ−ＭC -1ＢC ｜の値を用いて調節することができる。
【００２５】
処理装置は、不完全ブロック分解を段階的に行なう過程で、各レベルにおいて必要な逆行列ＭFF -1を、与えられた精度で近似的に求める。そして、その逆行列を利用して、次のレベル（粗いレベル）の分解で必要な行列ＡC を求め、ＡC のスパース性（sparsity）を、近似精度を考慮して決定する。また、不完全ブロック分解の各レベルで、行列ＡC の要素数がなるべく小さくなるように、言い換えれば、行列ＡFFの要素数がなるべく大きくなるように、集合Ｆの要素を決定する。
【００２６】
したがって、マルチレベル不完全ブロック分解では、図２に示すように、集合Ｆの決定、ＭFF -1の計算、およびＭC の非ゼロ要素の削減が、再帰的に繰り返される。（３）式において、ＡFFが対角行列であれば、その逆行列ＡFF -1は容易に求められる。この条件を少し緩めて、ＡFFが対角優位になるように集合Ｆを決定できれば、ＭFF -1は、反復計算により、任意の近似精度で求めることができる。ここで、ＡFFが対角優位であるとは、λを所定の定数として、ＡFFの要素が次の条件を満たすことを意味する。
【００２７】
【数３】
【００２８】
このとき、集合Ｆがなるべく大きくなるようにすれば、ＡCCが小さくなり、ＡC も小さくなるので、行列の分解が促進される。このような不完全ブロック分解により、Ｍ行列でない行列でも適切な近似精度で分解することができ、反復法の収束を速める前処理が可能になる。
【００２９】
図３および図４は、マルチレベル不完全ブロック分解を含む前処理のフローチャートである。再帰的アルゴリズムの各レベルにおいて、図３の処理を行うことで集合Ｆが決定され、粗いレベルのＡC を構成する変数番号が選び出される。
【００３０】
処理装置は、まず、対角優位なＡFFを抽出するために、分解の対象となる行列Ａの各行ベクトルを、最初に対角要素がくるように並べ替える（図３のステップＳ１１）。ただし、２番目以降のレベルでは、前のレベルの分解により得られたＡC が行列Ａとして用いられる。
【００３１】
連立一次方程式の係数行列Ａは、エルパック（Ｅｌｌｐａｃｋ）形式格納法により、処理装置のメモリに格納される。このエルパック形式格納法は、行ベクトルの非ゼロ要素を左方向に圧縮して格納し、非ゼロ要素があった列の番号を、整数型配列の対応する要素に格納する方法である。このとき、圧縮できなかった部分にはゼロを挿入し、整数型配列の対応する要素には、列番号の代わりに行番号を格納する。このように、スパース行列を圧縮して格納することで、メモリの使用量を節約することができる。
【００３２】
図５は、このようなエルパック形式格納法の例を示している。図５において、係数行列Ａの情報は、実数型２次元配列ＡＥと整数型２次元配列ＩＣＯＬに分けて格納される。配列ＡＥは、エルパック形式格納法で圧縮されたＡの行ベクトルを格納し、配列ＩＣＯＬは、Ａの非ゼロ要素の列番号を格納している。また、配列ＡＥのゼロ要素に対応する配列ＩＣＯＬの要素には、その要素の行番号が格納されている。
【００３３】
処理装置は、エルパック形式格納法で格納された行列の各行ベクトルの要素を並べ替えて、第１要素にＡの対角要素を格納し、他の要素については絶対値および対角要素からの距離を考慮して、それらの値が大きい順に格納する。また、これに対応してＩＣＯＬの要素も並べ替える。
【００３４】
ＩＣＯＬの（ｉ，ｊ）要素をｋi,j とすると、元の係数行列Ａの要素はａi,ki,jと表現できる。このとき、行ベクトルの要素の並べ替えは、以下の条件で行なわれる。
（１＋γ）｜ａi,ki,j+1｜≦｜ａi,ki,j｜ ｏｒ
｛（１−γ）｜ａi,ki,j+1｜≦｜ａi,ki,j｜ ａｎｄ
｜ｋi,j+1 −ｉ｜≦｜ｋi,j −ｉ｜｝
γは所定の定数であり、例えば、γ＝０．１５である。これにより、絶対値が同じ値の要素は、対角要素からの距離が大きいものほど左に位置するようになる。
【００３５】
次に、処理装置は、Ａの非ゼロ要素から構成されるグラフＧの頂点（ノード）の集合Ｖを生成する（ステップＳ１２）。このとき、ＡFFが対角優位を保つように、すべての頂点の集合から対角優位でない要素に対応する頂点を除いた集合をＶとし、除かれた頂点の集合を直接隣接集合Ｎの初期値とする。
【００３６】
まず、Ａの要素の添え字ｉを頂点（未知数）とするグラフを考え、頂点の対（ｉ，ｊ）のうち、ａi,j またはａj,i がゼロでないような（ｉ，ｊ）を抽出して、これらの頂点をエッジでリンクする。したがって、エッジ（ｉ，ｊ）は、Ａの非ゼロ要素に対応する。そして、グラフＧに属する頂点を識別するために、１次元配列にフラグを設定する。
【００３７】
次に、Ａから対角優位を崩すような要素を抽出し、直接隣接集合Ｎの初期値を決定する。このとき、並べ替えられたすべての行ベクトルについて、以下の判定を行毎に行なう。判定は、行ベクトルの２番目の要素から順に右の要素に関して行われ、すべての要素が終わるまで繰り返される。
そして、すべての頂点の集合から直接隣接集合Ｎに含まれる頂点を除いて、集合Ｖを生成する。
【００３８】
次に、処理装置は、ステップＳ１３〜Ｓ１７の処理により、集合Ｖから、集合Ｆとして用いられる最大独立集合を取り出す。ここでは、Ｖの部分集合Ｉのうち、Ｉの中のどんな２つの添え字を取り出しても、それらが上述のエッジ（ｉ，ｊ）に対応しないような部分集合（独立集合）を求める。そして、それらの独立集合のうち、要素数が最大なものを最大独立集合として選択する。この最大独立集合は、以下のアルゴリズムで求められる。
ＩＳＥＴ（Ｇ）：グラフＧから独立集合を取り出す手続き
ＳＥＧ（ｇｒａｐｈ（Ｇ），Ｖ）：グラフＧの部分集合のエッジの集合であって、各エッジの２つの頂点がＶに属するもの
Ｎ（Ｉ）：グラフＧのエッジの一方の頂点が独立集合Ｉに属するとき、Ｉに属さないもう一方の頂点の集合（直接隣接集合）
Ｅ：グラフＧのエッジの集合
Ｆ：最大独立集合（初期値はφ）
グラフＧは、エッジの総数をｍとして、２次元配列Ｇ（ｍ，２）で表現される。行列Ａの非ゼロ要素ａijに対して、頂点ｉ、ｊが共にグラフＧに属するとき、頂点番号ｉ、ｊが配列要素Ｇ（ｋ，１）、Ｇ（ｋ，２）にそれぞれ格納される（ｋ＝１，．．．，ｍ）。
【００３９】
ＩＳＥＴ（Ｇ）の処理では、例えば、乱数を利用して独立集合が生成される。この場合、各頂点番号に対して（０、１］の範囲の乱数を対応付け、配列要素Ｇ（ｋ，１）およびＧ（ｋ，２）に対応する２つの乱数を比較して、値が小さい方の頂点を独立集合に組み入れ、もう一方の頂点を独立集合から除外する。
【００４０】
このとき、各頂点に乱数を対応させるために、１次元配列に乱数を格納する。また、最大独立集合に属するか否かを示すフラグを格納する１次元配列（最大独立フラグ配列）を用意し、各頂点に対応するフラグの初期値として値“ｏｎ”を格納しておく。そして、乱数を比較した結果、大きい値に対応する頂点のフラグを値“ｏｆｆ”に変更する。こうして、フラグがｏｎになっている頂点の集合が、独立集合として取り出される。
【００４１】
以上のアルゴリズムを図３の処理ステップと対応付けて説明すると、次のようになる。まず、処理装置は、集合ＶよりグラフＧを生成し、各頂点に乱数を割り当てる（ステップＳ１３）。次に、２つの頂点の乱数の値を比較して、大きな値に対応する頂点を除き、残った頂点の集合を独立集合Ｉとする（ステップＳ１４）。
【００４２】
次に、集合Ｖ、Ｉ、およびＦを更新する（ステップＳ１５）。ここでは、集合Ｖから、独立集合Ｉおよび直接隣接集合Ｎ（Ｉ）を除いて、新たな集合Ｖを生成する。また、Ｎ（Ｉ）を直接隣接集合の和集合Ｎに加え、独立集合Ｉを最大独立集合Ｆに加える。直接隣接集合Ｎは、グラフのすべての頂点から最大独立集合Ｆの要素を除くことによっても、求めることができる。
【００４３】
次に、集合の更新を終了するか否かを判定する（ステップＳ１６）。ここでは、独立集合ＩがＮＵＬＬになるか、または、ＶまたはＥ（残りの頂点またはエッジ）が空になったときに、更新を終了するものとする。
【００４４】
更新を終了しない場合は、ステップＳ１３以降の処理を繰り返し、更新が終了すると、次に、集合Ｖを最大独立集合Ｆに加えて、最終的な最大独立集合Ｆを決定する（ステップＳ１７）。このとき、最大独立フラグ配列においてフラグがｏｎになっている頂点が、集合Ｆの変数番号を表している。また、直接隣接集合Ｎに属するためにグラフから除かれた頂点であることを示すフラグを格納する１次元配列（直接隣接フラグ配列）も用意され、対応する頂点のフラグがｏｎに設定される。
【００４５】
次に、処理装置は、得られた最大独立集合の変数番号の変換表を作成する（ステップＳ１８）。ここでは、グラフの頂点から最大独立集合に属する変数番号を識別するために、最大独立フラグ配列を参照し、最終的にフラグがｏｎの要素に、最大独立集合に属する頂点の通し番号を振り直す。そして、その番号を、最大独立フラグ配列のフラグの代わりに設定する。これにより、最大独立フラグ配列を利用して、元の変数番号の集合から最大独立集合の新たな変数番号への変換表が作成される。
【００４６】
次に、直接隣接集合の変数番号の変換表を作成する（ステップＳ１９）。この場合、最大独立集合と同様にして、直接隣接フラグ配列を参照し、フラグがｏｎの要素に、直接隣接集合に属する頂点の通し番号を振り直す。これにより、直接隣接フラグ配列を利用して、元の変数番号の集合から直接隣接集合の新たな変数番号への変換表が作成される。
【００４７】
次に、得られた変換表を利用して、ＡFFとＡCFを分離する（図４のステップＳ２０）。このとき、Ａの行ベクトルの要素に関して、変換表の新たな変数番号を参照し、最大独立集合に属する頂点だけを取り出して、それらに対応する行ベクトルの要素を左に詰めて配列ＡＥに格納する。また、これに応じて配列ＩＣＯＬを書き直す。そして、ＡＥとＩＣＯＬの行ベクトルをＡFFとＡCFに振り分けて、ＡFFとＡCFを分離する。
【００４８】
次に、ＡFFの近似逆行列ＦFFを反復法で求める（ステップＳ２１）。ＡFFは対角優位なので、その対角要素からなる対角行列ＤFFが存在する。したがって、ＤFFの逆行列ＤFF -1を初期解としてＭFF -1を反復法により求めれば、この計算は収束する。このとき、｜Ｉ−ＭFF -1ＡFF｜を反復的に計算して期待する近似精度でＭFF -1を求め、得られたＭFF -1をＦFFとする。具体的には、ＭFF,1 -1＝ＤFF -1とおいて、再帰的に以下の計算を行なう。
ＭFF,p -1＝ＭFF,p-1 -1＋（Ｉ−ＭFF,p-1 -1ＡFF）ＤFF -1 （６）
ρp ＝｜Ｉ−ＭFF,p -1ＡFF｜ （７）
そして、τF を所定の定数として、ρp ＜τF を満たしたところまでの近似を採用する。（６）式には、エルパック形式で格納されたスパース行列同士の行列積が現れるが、この計算は以下の手順で行われる。
【００４９】
一般に、Ａ＝Ａ＋αＢ×Ｃの形の行列演算では、まず、メモリの一時領域でｔ＝αＢ×Ｃを計算してから、Ａ＝Ａ＋ｔを計算する。ｔ＝αＢ×Ｃの計算では、エルパック形式で格納されたＢおよびＣの各々の行ベクトルの非ゼロ要素の最大数を、それぞれ、ｎｒｂおよびｎｒｃとして、一時的に第２次元の大きさがｎｒｂ×ｎｒｃの領域を確保する。次に、Ｃの列ベクトルに対してＢの各要素を掛ける。α＝１の場合の手順は、以下の通りである。
（１）ｔの配列ｔｍｐを計算する。
【００５０】
ｉｉ＝１
ｄｏ ｉ＝１，ｎｒｃ
ｄｏ ｋ＝１，ｎｒｂ
ｄｏ ｊ＝１，ｎ
ｔｍｐ（ｊ，ｉｉ）＝ｂ（ｊ，ｋ）＊ｃ（ｂｉｃｏｌ（ｊ，ｋ），ｉ）
ｔｍｐｉｃｏｌ（ｊ，ｉｉ）＝ｃｉｃｏｌ（ｂｉｃｏｌ（ｊ，ｋ），ｉ）
ｅｎｄｄｏ
ｉｉ＝ｉｉ＋１
ｅｎｄｄｏ
ｅｎｄｄｏ
（２）ｔｍｐの行ベクトルの要素を、それらが属する列ベクトルの番号順にソートする。
（３）Ａの行ベクトルも、同様にしてソートする。
（４）Ａ＝Ａ＋ｔｍｐを計算する。ここでは、行ベクトル毎に以下の処理を行う。まず、Ａ（ｉ，ｊ）のＩＣＯＬ（ｉ，ｊ）に一致するｔｍｐｉｃｏｌ（ｉ，ｋ）を、ｋを増加させながら探す。一致する要素が見つかれば、Ａ（ｉ，ｊ）にｔｍｐ（ｉ，ｋ）を加える。また、ｔｍｐｉｃｏｌ（ｉ，ｋ）と一致するＩＣＯＬ（ｉ，ｊ）がなければ、ｔｍｐ（ｉ，ｋ）およびｔｍｐｉｃｏｌ（ｉ，ｋ）を、ＡおよびＩＣＯＬの第ｉ行に新しい要素として追加する。
【００５１】
次に、前処理の計算で用いられるＡCFＦFFを計算し（ステップＳ２２）、ＡFCとＡCCを分離する（ステップＳ２３）。このとき、ステップＳ２０の処理と同様に、Ａの行ベクトルの要素に関して、変換表の新たな変数番号を参照し、直接隣接集合に属する頂点だけを取り出し、配列ＩＣＯＬを書き直す。そして、行ベクトルをＡFCとＡCCに振り分ける。
【００５２】
次に、ＡC ＝ＡCC−ＡCFＦFFＡFCを計算する（ステップＳ２４）。ここでは、ステップＳ２１の行列積の計算と同様にして、ＡCFＦFFＡFCをｔｍｐに格納し、ＡCFＦFFＡFCをＡCCから差し引く。このとき、ＡC の各要素の絶対値を判定しながら、与えられた値より小さな要素を落として、ＭC ＝ＡC −ＲC によりＡC を近似し、ＡC のスパース性を保つ。そして、ＭC をＡC の計算結果として採用する。
【００５３】
このとき、例えば、τを所定の定数として、（ＡC の要素の絶対値）＜τ×（要素の属する行ベクトルの各要素の絶対値の和）という条件を満たすか否かを判定し、この条件を満たす要素を落とすことにする。具体的には、ＲC を決定するための条件は、｜ＭC -1ＲC ｜1 ＜τで与えられる。ただし、ｎ×ｎの行列Ｂ＝（ｂi,j ）に対して、｜Ｂ｜1 は次式で定義される。
【００５４】
【数４】
【００５５】
ここで、未知行列ＭC を以下の対角行列で近似して考えることにする。
【００５６】
【数５】
【００５７】
ここで、ｓｉｇｎ（ｘ）は、ｘ≧０のとき１となり、ｘ＜０のとき−１となる。このような近似によれば、ＡC の要素ａi,j が次式を満たすときに、その要素を落とすことになる。
｜ａi,j ｜＜τ×ｍi,i （１０）
実際には、ＡC の行ベクトルの各要素の絶対値の和を求めて、その行の各要素に関して、左から順に（１０）式の条件を判定し、条件を満たせばその要素を削除する。そして、次の要素を詰めて格納する。
【００５８】
次に、得られたＡC が分解の最も粗いレベル（最終レベル）に対応するか否かを判定する（ステップＳ２５）。最終レベルに対応しなければ、得られたＡC をＡとして用いて、図３のステップＳ１１以降の処理を繰り返す。そして、最終レベルに達したら、行列の分解を終了し、その結果を利用して前処理の計算を行い（ステップＳ２６）、処理を終了する。前処理の計算においては、上述したように、各レベルの不完全ブロック分解を利用して、行列ベクトル積Ｍ-1ｑが求められる。
【００５９】
次に、図６から図９までを参照しながら、上述した前処理を並列に行なう構成について説明する。
図６は、対称型マルチプロセッサ（Symmetric Multi-Processors，ＳＭＰ）システム（並列計算機）の構成図である。図６のマルチプロセッサシステムは、複数のプロセッサ１１、入力装置１２、出力装置１３、および共有メモリ１４を備え、それらはシステムバス１５により接続されている。入力装置１２は、ユーザからの指示や情報の入力に用いられ、出力装置１３は、ユーザへの問い合わせや処理結果の出力に用いられる。
【００６０】
プロセッサ１１は、キャッシュ付きのＣＰＵ（中央処理装置）に対応し、共有メモリ１４を介して互いに結合されている。この場合、並列処理のアプリケーションプログラムは、スレッド（thread）と呼ばれる処理単位を利用して、ＳＭＰの並列拡張言語であるＯｐｅｎＭＰにより記述される。
【００６１】
ＯｐｅｎＭＰでは、変数にＳＨＡＲＥＤおよびＰＲＩＶＡＴＥの属性を付与することができる。このうち、ＳＨＡＲＥＤ属性を持つ変数は、スレッド間で共通に参照・書き込みを行うことができる。また、ＯｐｅｎＭＰでは、ＢＡＲＲＩＥＲ同期を利用して、各スレッドの番号を判定しながら、並列処理を制御することができる。
【００６２】
このシステムでは、処理対象の頂点番号がプロセッサの台数で分割される。そして、図７に示すように、配列ＡＥおよびＩＣＯＬにエルパック形式で格納された行列Ａにおいて、各プロセッサに割り当てられた番号に関する要素の計算を、各スレッドが担当する。ここでは、互いに異なるプロセッサ上に生成されたスレッドＴ１、Ｔ２、Ｔ３、およびＴ４により、行列Ａが分割されている。
【００６３】
ＡFF、ＡCF、ＡFC、およびＡCCの各スレッドの担当部分は、スレッド毎に連続な共用領域を作業領域として用いて格納される。このとき、図８に示すように、行列の非ゼロ要素が存在する部分のバンド幅（上バンド幅および下バンド幅）が求められ、各スレッドの担当部分のバンド幅が決められる。
【００６４】
また、図３のステップＳ１３で各頂点に割り当てられる乱数を格納する配列、最大独立フラグ配列、直接隣接フラグ配列、および変数番号の変換表等の格納領域は、図９に示すように設定される。図９では、各スレッドが担当する部分の格納領域２２の左側に、下バンド幅分の余分な領域（オーバラップ領域）２１が設けられ、右側に上バンド幅分のオーバラップ領域２３が設けられている。
【００６５】
ステップＳ１３〜Ｓ１７の処理においては、最大独立集合を決定するアルゴリズムを以下の方法で並列化する。ここでは、フラグの“ｏｎ”と“ｏｆｆ”をそれぞれ整数の“１”と“０”で表し、組み込み関数ＭＩＮを利用して設定する。最大独立フラグ配列の初期値としては、直接隣接集合の頂点以外をｏｎとし、直接隣接集合として除かれた頂点をｏｆｆとする。
【００６６】
まず、各スレッドは、乱数を生成して、グラフを構成する頂点のオーバラップ領域以外の領域（各スレッドが担当する領域）に格納する。その後、オーバラップ領域に、その部分を担当する他のスレッドの領域から乱数をコピーして、値を設定し直す。これにより、各スレッドで用いる頂点の乱数が、連続領域でローカルに保持される。
【００６７】
次に、各スレッドは、担当するグラフを表現するエッジの頂点対をすべて処理する。処理対象の頂点は、各スレッドが担当する部分の頂点番号に上バンド幅分および下バンド幅分の頂点番号を加えた範囲に存在する。そして、頂点対に対応する２つの乱数を比較し、最大独立フラグ配列において、値が大きい方の頂点のフラグをｏｆｆにする。
【００６８】
各スレッドでの設定が終了すると、各スレッドは、オーバラップ領域の頂点を担当する他のスレッドからフラグ情報をコピーして、最大独立フラグ配列のオーバラップ領域に格納する。このとき、関数ＭＩＮを用いてスレッド間で整合性をとり、他のスレッドでｏｆｆになっている部分はｏｆｆにする。
【００６９】
次に、各スレッドは、もう一度グラフのエッジを走査し、一方の頂点が最大独立集合に属する場合、直接隣接フラグ配列において、もう一方の頂点のフラグをｏｎにする。また、最大独立フラグ配列と同様にして、スレッド間で整合性をとるように、直接隣接フラグ配列のオーバラップ領域の値を設定する。
【００７０】
次に、各スレッドは、グラフを再度走査し、最大独立集合と直接隣接集合に属さないエッジのみを、グラフを構成するエッジとして残す。このとき、グラフを構成するエッジだけを詰め、最大独立集合の要素数をカウントして、それがゼロか否かをチェックする。また、頂点およびエッジの数をカウントして、それらがゼロか否かをチェックする。最大独立集合の要素数、頂点およびエッジの数が共にゼロでなければ、新たなグラフに関して、乱数の生成以降の処理を繰り返す。
【００７１】
そして、最大独立集合の要素数またはエッジの数がゼロになると、各スレッドは、繰り返しを終了し、直接隣接フラグ配列において、エッジを構成する頂点のフラグをｏｎにする。そして、スレッド間で格納領域の整合性をとる。
【００７２】
次に、ステップＳ１８において、アプリケーションは、各スレッドが担当する頂点（オーバラップ領域は除く）のうち、最大独立集合に属するものの数をカウントする。そして、スレッドの総数の大きさの１次元配列（共有領域）に各スレッドのカウント値を格納し、この情報を用いて、各スレッドの最大独立集合の頂点に、全体での通し番号を振る。このとき、通し番号は、正整数で１から順に振られる。
【００７３】
次に、各スレッドは、最大独立フラグ配列のｏｎの代わりに、その頂点に対応する通し番号を書き込み、変数番号の変換表を作成する。この変換表において、正整数に対応する変数番号は、最大独立集合に属する頂点に対応し、０に対応する変数番号は、直接隣接集合に属する頂点に対応する。そして、オーバラップ領域に対応する他のスレッドから、変換表の一部をコピーする。ステップＳ１９においても、同様にして、直接隣接集合の変数番号の変換表が作成される。
【００７４】
次に、図４のステップＳ２０において、各スレッドは、並列に作業領域を確保し、最大独立集合の変換表を用いてＡFFとＡCFを分離する。この場合、Ａの要素ａi,j のうち、ｊが最大独立集合に属するもののみを、ＡＥ内で左に詰めて格納し、ＩＣＯＬの格納値をその頂点の新たな変数番号に書き換える。
【００７５】
次に、各行ベクトルの非ゼロ要素の個数の最大値Ｎｒを求め、各スレッドが担当する最大独立集合および直接隣接集合の要素数とＮｒに基づいて、各スレッドの作業領域にＡFFとＡCFの格納領域（ＡＥおよびＩＣＯＬ）を確保する。そして、ＡのＡＥおよびＩＣＯＬの行ベクトルをＡFFとＡCFの格納領域に振り分ける。ステップＳ２３においても、同様にして、直接隣接集合の変換表を用いてＡFCとＡCCが分離される。
【００７６】
また、ステップＳ２１およびＳ２４においては、スパース行列積の計算が並列に行われる。各スレッドは、図７および図８と同様に、処理対象の各行列の対応部分を保持している。
【００７７】
上述したｔｍｐの計算において、最も内側のｄｏループの処理は、各スレッドに割り当てられた部分に対してのみ行われる。また、必要なＣの列ベクトルに関しても、上バンド幅と下バンド幅のオーバラップ領域を有する作業領域を用意し、各スレッドが保持する部分と他のスレッドが保持する部分をその作業領域にコピーしてから、各スレッドで並列に計算する。
【００７８】
また、ステップＳ２６のｐ＝Ｍ-1ｑの並列計算において、アプリケーションは、ベクトルｑF およびｑC の各々に対して、各スレッドが担当する部分の格納領域と上下バンド幅のオーバラップ領域を用意し、変換表を利用して、ベクトルｑからｑF およびｑC を分離する。そして、各スレッドは、割り当てられた行列部分を用いて、行列ベクトル積を並列に計算する。このとき、他のスレッドの処理結果を参照するために、オーバラップ領域を用いる。
【００７９】
ところで、本実施形態の処理装置は、並列計算機に限らず、図１０に示すような任意の情報処理装置（コンピュータ）を用いて構成することもできる。図１０の情報処理装置は、ＣＰＵ（中央処理装置）３１、メモリ３２、入力装置３３、出力装置３４、外部記憶装置３５、媒体駆動装置３６、およびネットワーク接続装置３７を備え、それらはバス３８により互いに接続されている。
【００８０】
メモリ３２は、例えば、ＲＯＭ（read only memory）、ＲＡＭ（random access memory）等を含み、処理に用いられるプログラムとデータを格納する。ＣＰＵ３１は、メモリ３２を利用してプログラムを実行することにより、必要な処理を行う。上述した連立１次方程式の反復法および前処理のプログラムは、インストラクションの集合として、メモリ３２の特定のプログラムコードセグメントに格納される。
【００８１】
入力装置３３は、例えば、キーボード、ポインティングデバイス、タッチパネル等であり、ユーザからの指示や情報の入力に用いられる。出力装置３４は、例えば、ディスプレイ、プリンタ、スピーカ等であり、ユーザへの問い合わせや処理結果の出力に用いられる。
【００８２】
外部記憶装置３５は、例えば、磁気ディスク装置、光ディスク装置、光磁気ディスク（magneto-optical disk）装置、テープ装置等である。情報処理装置は、この外部記憶装置３５に、上述のプログラムとデータを保存しておき、必要に応じて、それらをメモリ３２にロードして使用する。
【００８３】
媒体駆動装置３６は、可搬記録媒体３９を駆動し、その記録内容にアクセスする。可搬記録媒体３９としては、メモリカード、フロッピーディスク、ＣＤ−ＲＯＭ（compact disk read only memory ）、光ディスク、光磁気ディスク等、任意のコンピュータ読み取り可能な記録媒体が用いられる。ユーザは、この可搬記録媒体３９に上述のプログラムとデータを格納しておき、必要に応じて、それらをメモリ３２にロードして使用する。
【００８４】
ネットワーク接続装置３７は、ＬＡＮ（local area network）等の任意の通信ネットワークに接続され、通信に伴うデータ変換を行う。また、情報処理装置は、上述のプログラムとデータをネットワーク接続装置３７を介して、サーバ等の他の装置から受け取り、必要に応じて、それらをメモリ３２にロードして使用する。
【００８５】
図１１は、図１０の情報処理装置にプログラムとデータを供給することのできるコンピュータ読み取り可能な記録媒体を示している。可搬記録媒体３９やサーバ４０のデータベース４１に保存されたプログラムとデータは、メモリ３２にロードされる。そして、ＣＰＵ３１は、そのデータを用いてそのプログラムを実行し、必要な処理を行う。このとき、サーバ４０は、プログラムとデータを伝送する伝搬信号を生成し、ネットワーク上の任意の伝送媒体を介して、情報処理装置に送信する。
（付記１） 連立一次方程式を解く反復法のための前処理を、マルチレベル不完全ブロック分解により行なう処理装置であって、
前記不完全ブロック分解のあるレベルにおける係数行列を格納する係数格納手段と、
前記係数格納手段に格納された係数行列の要素のうち、取り除かれる変数の変数番号を行番号および列番号とする要素からなるブロック行列が対角優位になるように、該取り除かれる変数の変数番号の集合を決定する決定手段と、
決定された変数番号の集合に関する情報を格納する集合格納手段と、
前記集合格納手段に格納された情報に基づいて、前記ブロック行列の近似逆行列を求める逆行列手段と、
前記係数格納手段に格納された係数行列の要素のうち、残される変数の変数番号を行番号および列番号とする要素からなるブロック行列から、粗いレベルの係数行列を求める係数行列手段と、
次の不完全ブロック分解を行うために、前記粗いレベルの係数行列を前記係数格納手段に格納する分解手段と、
前記近似逆行列を用いて、前記反復法における行列ベクトル積を計算する計算手段と
を備えることを特徴とする処理装置。
（付記２） 前記決定手段は、前記係数格納手段に格納された係数行列の要素の添え字の集合の部分集合である直接隣接集合の初期値を設定し、該直接隣接集合の初期値と乱数を用いて添え字の最大独立集合を求めることで、前記変数番号の集合を決定することを特徴とする付記１記載の処理装置。
（付記３） 前記逆行列手段は、与えられた近似精度で反復計算により前記近似逆行列を求めることを特徴とする付記１記載の処理装置。
（付記４） 前記係数行列手段は、前記粗いレベルの係数行列の非ゼロ要素を削減して、該粗いレベルの係数行列のスパース性を保持することを特徴とする付記１記載の処理装置。
（付記５） 複数のプロセッサを有し、連立一次方程式を解く反復法のための前処理を、マルチレベル不完全ブロック分解により行なう並列計算機であって、
前記複数のプロセッサが前記不完全ブロック分解のあるレベルにおける係数行列を並列に処理できるように、該係数行列を分割して格納する係数格納手段と、前記係数格納手段に格納された係数行列の要素のうち、取り除かれる変数の変数番号を行番号および列番号とする要素からなるブロック行列が対角優位になるように、該取り除かれる変数の変数番号の集合を並列に決定する決定手段と、
前記複数のプロセッサが決定された変数番号の集合に関する情報を並列に処理できるように、該変数番号の集合に関する情報を分割して格納する集合格納手段と、
前記集合格納手段に格納された情報に基づいて、前記ブロック行列の近似逆行列を並列に求める逆行列手段と、
前記係数格納手段に格納された係数行列の要素のうち、残される変数の変数番号を行番号および列番号とする要素からなるブロック行列から、粗いレベルの係数行列を並列に求める係数行列手段と、
次の不完全ブロック分解を行うために、前記粗いレベルの係数行列を前記係数格納手段に分割して格納する分解手段と、
前記近似逆行列を用いて、前記反復法における行列ベクトル積を並列に計算する計算手段と
を備えることを特徴とする並列計算機。
（付記６） 連立一次方程式を解く反復法のための前処理を、マルチレベル不完全ブロック分解により行なうコンピュータのためのプログラムを記録した記録媒体であって、該プログラムは、
前記不完全ブロック分解のあるレベルにおける係数行列を格納し、
格納された係数行列の要素のうち、取り除かれる変数の変数番号を行番号および列番号とする要素からなるブロック行列が対角優位になるように、該取り除かれる変数の変数番号の集合を決定し、
決定された変数番号の集合に関する情報を格納し、
格納された情報に基づいて、前記ブロック行列の近似逆行列を求め、
前記格納された係数行列の要素のうち、残される変数の変数番号を行番号および列番号とする要素からなるブロック行列から、粗いレベルの係数行列を求め、次の不完全ブロック分解を行うために、前記粗いレベルの係数行列を前記あるレベルにおける係数行列として格納し、
前記近似逆行列を用いて、前記反復法における行列ベクトル積を計算する
処理を前記コンピュータに実行させることを特徴とするコンピュータ読み取り可能な記録媒体。
（付記７） 連立一次方程式を解く反復法のための前処理を、マルチレベル不完全ブロック分解により行なうコンピュータにプログラムを伝送する伝搬信号であって、該プログラムは、
前記不完全ブロック分解のあるレベルにおける係数行列を格納し、
格納された係数行列の要素のうち、取り除かれる変数の変数番号を行番号および列番号とする要素からなるブロック行列が対角優位になるように、該取り除かれる変数の変数番号の集合を決定し、
決定された変数番号の集合に関する情報を格納し、
格納された情報に基づいて、前記ブロック行列の近似逆行列を求め、
前記格納された係数行列の要素のうち、残される変数の変数番号を行番号および列番号とする要素からなるブロック行列から、粗いレベルの係数行列を求め、
次の不完全ブロック分解を行うために、前記粗いレベルの係数行列を前記あるレベルにおける係数行列として格納し、
前記近似逆行列を用いて、前記反復法における行列ベクトル積を計算する
処理を前記コンピュータに実行させることを特徴とする伝搬信号。
【００８６】
【発明の効果】
本発明によれば、係数行列がＭ行列でないような連立一次方程式の反復法において、計算の収束を速めるような前処理を行うことが可能になる。したがって、従来の前処理では収束しなかったような問題を、効率良く解くことができる。
【図面の簡単な説明】
【図１】本発明の処理装置の原理図である。
【図２】マルチレベル不完全ブロック分解を示す図である。
【図３】前処理のフローチャート（その１）である。
【図４】前処理のフローチャート（その２）である。
【図５】エルパック形式格納法を示す図である。
【図６】マルチプロセッサシステムの構成図である。
【図７】配列の分割を示す図である。
【図８】バンド幅を示す図である。
【図９】格納領域を示す図である。
【図１０】情報処理装置の構成図である。
【図１１】記録媒体を示す図である。
【符号の説明】
１ 係数格納手段
２ 決定手段
３ 集合格納手段
４ 逆行列手段
５ 係数行列手段
６ 分解手段
７ 計算手段
１１ プロセッサ
１２、３３ 入力装置
１３、３４ 出力装置
１４ 共有メモリ
１５ システムバス
２１、２３ オーバラップ領域
２２ 格納領域
３１ ＣＰＵ
３２ メモリ
３５ 外部記憶装置
３６ 媒体駆動装置
３７ ネットワーク接続装置
３８ バス
３９ 可搬記録媒体
４０ サーバ
４１ データベース
【発明の属する技術分野】
本発明は、連立一次方程式の解を求める反復法に係り、さらに詳しくは、反復法における前処理（preconditioning ）を行う処理装置に関する。
【０００２】
【従来の技術】
物理現象の解析においては、しばしば、大規模な行列を係数行列とする次のような連立一次方程式を解く必要が生じる。
Ａｘ＝ｂ （１）
ここで、Ａはｎ×ｎの係数行列、ｘはｎ次元の変数ベクトル、ｂはｎ次元の定数ベクトルである。ｎの値が１０6 程度になることも珍しくない。
【０００３】
このような大規模な連立一次方程式は、気象予測、原子炉設計、半導体の回路解析、航空工学における流体解析、構造物の構造解析等の多くの科学技術計算に用いられる。したがって、大規模な連立一次方程式を効率よく高速に解くことは、科学技術計算の重要な問題の１つである。
【０００４】
コンピュータを用いて（１）式を解く１つの方法として、ＡをＬＵ分解するガウスの消去法に基づいた直接法がある。しかし、Ａが大きなスパース行列の場合、非零要素が各行に数個しかないこともあり、計算コストや記憶領域の面で無駄が多い。そこで、このような場合には、単純な行列ベクトル積を繰り返して近似解を求める反復法が多く用いられている。
【０００５】
スパース行列の反復法の収束を促進する方法として、前処理が用いられる。この処理では、適当な前処理行列（preconditioner）Ｍ-1を用いて、（１）式が次のように変更される。
Ｍ-1Ａｘ＝Ｍ-1ｂ （２）
そして、（２）式に対して反復法を適用することで、（１）式の近似解が求められる。このとき、ＭはＡの不完全ブロック分解（block incomplete factorization，ＢＩＦ）に対応し、Ｍ-1はＭの逆行列を表す。このように、Ａの不完全ブロック分解を利用して前処理を行う１つの方法として、代数的マルチレベル反復法（algebraic multilevel iteration method ）が知られている。
【０００６】
【発明が解決しようとする課題】
しかしながら、上述した従来の代数的マルチレベル反復法による前処理には、次のような問題がある。
【０００７】
一般に、代数的マルチレベル反復法による前処理は、係数行列ＡがＭ行列（M-matrix）であるような問題に対して適用でき、反復法の収束を速める効果があることが知られている。ここで、ＡがＭ行列であるとは、Ａ＝（ａij）が正則行列であり、さらに次の２つの条件を満たすことを意味する。
（１）ａij≦０（ｉ≠ｊ）（Ａの非対角要素が０もしくは負）
（２）Ａ-1≧０（Ａの逆行列Ａ-1の全要素が非負）
しかし、ＡがＭ行列でない場合は、代数的マルチレベル反復法による前処理はあまり効果がなく、計算が収束しないという問題がある。
【０００８】
本発明の課題は、係数行列がＭ行列でないような連立一次方程式の反復法において、計算の収束を速めるような前処理を行う処理装置を提供することである。
【０００９】
図１は、本発明の処理装置の原理図である。図１の処理装置は、係数格納手段１、決定手段２、集合格納手段３、逆行列手段４、係数行列手段５、分解手段６、および計算手段７を備え、係数行列がＭ行列でない連立一次方程式を解く反復法のための前処理を、マルチレベル不完全ブロック分解により行なう。
【００１０】
係数格納手段１は、不完全ブロック分解のあるレベルにおける係数行列を格納する。決定手段２は、係数格納手段１に格納された係数行列の要素のうち、取り除かれる変数の変数番号を行番号および列番号とする要素からなるブロック行列が対角優位（diagonal dominant ）になるように、それらの取り除かれる変数の変数番号の集合を決定する。このとき、決定手段２は、係数格納手段１に格納された係数行列の要素のうち、対角優位な各要素の添え字に乱数を割り当て、２つの添え字に割り当てられた乱数を比較して大きい方に対応する添え字を取り除く処理を繰り返すことで、残された添え字からなる最大独立集合を求め、その最大独立集合を変数番号の集合とする。そして、集合格納手段３は、決定された変数番号の集合に関する情報を格納する。
【００１１】
また、逆行列手段４は、集合格納手段３に格納された情報に基づいて、上記ブロック行列の近似逆行列を求める。係数行列手段５は、係数格納手段１に格納された係数行列の要素のうち、残される変数の変数番号を行番号および列番号とする要素からなるブロック行列から、次のレベルの係数行列を求める。そして、分解手段６は、次の不完全ブロック分解を行うために、次のレベルの係数行列を係数格納手段１に格納する。
【００１２】
また、逆行列手段４は、集合格納手段３に格納された情報に基づいて、上記ブロック行列の近似逆行列を求める。係数行列手段５は、係数格納手段１に格納された係数行列の要素のうち、残される変数の変数番号を行番号および列番号とする要素からなるブロック行列から、次のレベルの係数行列を求める。そして、分解手段６は、次の不完全ブロック分解を行うために、次のレベルの係数行列を係数格納手段１に格納する。
【００１３】
次に、逆行列手段４は、その集合の情報に基づいて、取り除かれる変数の係数からなるブロック行列の逆行列を近似的に求め、係数行列手段５は、上記係数行列の情報を参照して、処理対象として残される変数の係数からなるブロック行列から、次のレベルの係数行列を計算する。そして、分解手段６が、得られた係数行列を分解される係数行列として係数格納手段１に格納することで、次のレベルの不完全ブロック分解が再帰的に行われる。
【００１４】
計算手段７は、再帰的な不完全ブロック分解により得られた各レベルの近似逆行列を用いて、反復法で必要となる行列ベクトル積を再帰的に計算することで、反復法のための前処理を行う。
【００１５】
対角優位な行列の場合、後述するように、与えられた近似精度で反復計算により、その逆行列を求めることができる。したがって、取り除かれる変数の係数からなるブロック行列が対角優位になるように、それらの変数の変数番号の集合を決定することで、任意のスパース行列を適切な近似精度で分解することができる。このような不完全ブロック分解によれば、Ｍ行列でない係数行列に対しても、反復法の収束を速める前処理が可能になる。
【００１６】
例えば、図１の係数格納手段１および集合格納手段３は、後述する図６の共有メモリ１４または図１０のメモリ３２に対応し、図１の決定手段２、逆行列手段４、係数行列手段５、分解手段６、および計算手段７は、図６のプロセッサ１１または図１０のＣＰＵ（中央処理装置）３１に対応する。
【００１７】
また、取り除かれる変数の変数番号の集合は、後述する集合Ｆに対応し、それらの変数番号を行番号および列番号とする要素からなるブロック行列は、後述する（３）式の行列ＡFFに対応する。また、近似逆行列は、後述する（４）式のＭFF -1または図４のＦFFに対応し、残される変数の変数番号を行番号および列番号とする要素からなるブロック行列は、（３）式の行列ＡCCに対応し、次のレベルの係数行列は、図４の行列ＡC または（４）式のＭC に対応する。
【００１８】
【発明の実施の形態】
以下、図面を参照しながら、本発明の実施の形態を詳細に説明する。
一般に、前処理を施した反復法においては、（２）式のＭ-1とｎ次元ベクトルｑの行列ベクトル積Ｍ-1ｑを繰り返し計算する必要が生じる。本実施形態の処理装置は、代数的マルチレベル反復法による前処理において、連立一次方程式の係数行列のマルチレベル不完全ブロック分解を行い、各レベルで得られた結果を利用してＭ-1ｑを計算する。
【００１９】
まず、連立方程式の変数番号の集合｛１，．．．，ｎ｝から適当に選択された変数番号の集合をＣとし、その補集合をＦとすると、ｎ×ｎの行列Ａは以下の形に変換できる。
【００２０】
【数１】
ここで、集合Ｃは、残される変数を表す変数番号の集合に対応し、集合Ｆは、取り除かれる変数を表す変数番号の集合に対応する。また、ブロック行列ＡCCの要素の行番号および列番号は、集合Ｃの要素により表され、ＡCCの行および列の数は、集合Ｃの要素の数に一致する。同様に、ブロック行列ＡFFの要素の行番号および列番号は、集合Ｆの要素により表され、ＡFFの行および列の数は、集合Ｆの要素の数に一致する。
【００２２】
また、ＡFFの近似逆行列をＭFF -1とし、ＡC ：＝ＡCC−ＡCFＭFF -1ＡFCとすると、ＡC は、例えば、ＢC ＝ＡC −ＲC により近似することができる。ただし、ＲC は、近似により落とされる非ゼロ要素を含む行列である。ＡC の近似をＭC と書くことにすると、Ａの不完全ブロック分解Ｍは以下のように書ける。
【００２３】
【数２】
そして、ＭC の不完全ブロック分解を、粗いレベルからさらに粗いレベルへと再帰的に繰り返すことにより、各レベルにおいて、（４）式のような不完全ブロック分解を得ることができる。また、ｐ：＝Ｍ-1ｑ＝ＢＩＦ（０，ｑ）は、（４）式のブロック行列を用いて再帰的に計算することができる。最初のレベルを０とし、最も粗いレベルをＬとすると、この行列ベクトル積の計算アルゴリズムは、以下のようになる。
【００２５】
処理装置は、不完全ブロック分解を段階的に行なう過程で、各レベルにおいて必要な逆行列ＭFF -1を、与えられた精度で近似的に求める。そして、その逆行列を利用して、次のレベル（粗いレベル）の分解で必要な行列ＡC を求め、ＡC のスパース性（sparsity）を、近似精度を考慮して決定する。また、不完全ブロック分解の各レベルで、行列ＡC の要素数がなるべく小さくなるように、言い換えれば、行列ＡFFの要素数がなるべく大きくなるように、集合Ｆの要素を決定する。
【００２６】
したがって、マルチレベル不完全ブロック分解では、図２に示すように、集合Ｆの決定、ＭFF -1の計算、およびＭC の非ゼロ要素の削減が、再帰的に繰り返される。（３）式において、ＡFFが対角行列であれば、その逆行列ＡFF -1は容易に求められる。この条件を少し緩めて、ＡFFが対角優位になるように集合Ｆを決定できれば、ＭFF -1は、反復計算により、任意の近似精度で求めることができる。ここで、ＡFFが対角優位であるとは、λを所定の定数として、ＡFFの要素が次の条件を満たすことを意味する。
【００２７】
【数３】
このとき、集合Ｆがなるべく大きくなるようにすれば、ＡCCが小さくなり、ＡC も小さくなるので、行列の分解が促進される。このような不完全ブロック分解により、Ｍ行列でない行列でも適切な近似精度で分解することができ、反復法の収束を速める前処理が可能になる。
【００２９】
図３および図４は、マルチレベル不完全ブロック分解を含む前処理のフローチャートである。再帰的アルゴリズムの各レベルにおいて、図３の処理を行うことで集合Ｆが決定され、粗いレベルのＡC を構成する変数番号が選び出される。
【００３０】
処理装置は、まず、対角優位なＡFFを抽出するために、分解の対象となる行列Ａの各行ベクトルを、最初に対角要素がくるように並べ替える（図３のステップＳ１１）。ただし、２番目以降のレベルでは、前のレベルの分解により得られたＡC が行列Ａとして用いられる。
【００３１】
連立一次方程式の係数行列Ａは、エルパック（Ｅｌｌｐａｃｋ）形式格納法により、処理装置のメモリに格納される。このエルパック形式格納法は、行ベクトルの非ゼロ要素を左方向に圧縮して格納し、非ゼロ要素があった列の番号を、整数型配列の対応する要素に格納する方法である。このとき、圧縮できなかった部分にはゼロを挿入し、整数型配列の対応する要素には、列番号の代わりに行番号を格納する。このように、スパース行列を圧縮して格納することで、メモリの使用量を節約することができる。
【００３２】
図５は、このようなエルパック形式格納法の例を示している。図５において、係数行列Ａの情報は、実数型２次元配列ＡＥと整数型２次元配列ＩＣＯＬに分けて格納される。配列ＡＥは、エルパック形式格納法で圧縮されたＡの行ベクトルを格納し、配列ＩＣＯＬは、Ａの非ゼロ要素の列番号を格納している。また、配列ＡＥのゼロ要素に対応する配列ＩＣＯＬの要素には、その要素の行番号が格納されている。
【００３３】
処理装置は、エルパック形式格納法で格納された行列の各行ベクトルの要素を並べ替えて、第１要素にＡの対角要素を格納し、他の要素については絶対値および対角要素からの距離を考慮して、それらの値が大きい順に格納する。また、これに対応してＩＣＯＬの要素も並べ替える。
【００３４】
ＩＣＯＬの（ｉ，ｊ）要素をｋi,j とすると、元の係数行列Ａの要素はａi,ki,jと表現できる。このとき、行ベクトルの要素の並べ替えは、以下の条件で行なわれる。
（１＋γ）｜ａi,ki,j+1｜≦｜ａi,ki,j｜ ｏｒ
｛（１−γ）｜ａi,ki,j+1｜≦｜ａi,ki,j｜ ａｎｄ
｜ｋi,j+1 −ｉ｜≦｜ｋi,j −ｉ｜｝
γは所定の定数であり、例えば、γ＝０．１５である。これにより、絶対値が同じ値の要素は、対角要素からの距離が大きいものほど左に位置するようになる。
【００３５】
次に、処理装置は、Ａの非ゼロ要素から構成されるグラフＧの頂点（ノード）の集合Ｖを生成する（ステップＳ１２）。このとき、ＡFFが対角優位を保つように、すべての頂点の集合から対角優位でない要素に対応する頂点を除いた集合をＶとし、除かれた頂点の集合を直接隣接集合Ｎの初期値とする。
【００３６】
まず、Ａの要素の添え字ｉを頂点（未知数）とするグラフを考え、頂点の対（ｉ，ｊ）のうち、ａi,j またはａj,i がゼロでないような（ｉ，ｊ）を抽出して、これらの頂点をエッジでリンクする。したがって、エッジ（ｉ，ｊ）は、Ａの非ゼロ要素に対応する。そして、グラフＧに属する頂点を識別するために、１次元配列にフラグを設定する。
【００３７】
次に、Ａから対角優位を崩すような要素を抽出し、直接隣接集合Ｎの初期値を決定する。このとき、並べ替えられたすべての行ベクトルについて、以下の判定を行毎に行なう。判定は、行ベクトルの２番目の要素から順に右の要素に関して行われ、すべての要素が終わるまで繰り返される。
【００３８】
次に、処理装置は、ステップＳ１３〜Ｓ１７の処理により、集合Ｖから、集合Ｆとして用いられる最大独立集合を取り出す。ここでは、Ｖの部分集合Ｉのうち、Ｉの中のどんな２つの添え字を取り出しても、それらが上述のエッジ（ｉ，ｊ）に対応しないような部分集合（独立集合）を求める。そして、それらの独立集合のうち、要素数が最大なものを最大独立集合として選択する。この最大独立集合は、以下のアルゴリズムで求められる。
ＩＳＥＴ（Ｇ）：グラフＧから独立集合を取り出す手続き
ＳＥＧ（ｇｒａｐｈ（Ｇ），Ｖ）：グラフＧの部分集合のエッジの集合であって、各エッジの２つの頂点がＶに属するもの
Ｎ（Ｉ）：グラフＧのエッジの一方の頂点が独立集合Ｉに属するとき、Ｉに属さないもう一方の頂点の集合（直接隣接集合）
Ｅ：グラフＧのエッジの集合
Ｆ：最大独立集合（初期値はφ）
【００３９】
ＩＳＥＴ（Ｇ）の処理では、例えば、乱数を利用して独立集合が生成される。この場合、各頂点番号に対して（０、１］の範囲の乱数を対応付け、配列要素Ｇ（ｋ，１）およびＧ（ｋ，２）に対応する２つの乱数を比較して、値が小さい方の頂点を独立集合に組み入れ、もう一方の頂点を独立集合から除外する。
【００４０】
このとき、各頂点に乱数を対応させるために、１次元配列に乱数を格納する。また、最大独立集合に属するか否かを示すフラグを格納する１次元配列（最大独立フラグ配列）を用意し、各頂点に対応するフラグの初期値として値“ｏｎ”を格納しておく。そして、乱数を比較した結果、大きい値に対応する頂点のフラグを値“ｏｆｆ”に変更する。こうして、フラグがｏｎになっている頂点の集合が、独立集合として取り出される。
【００４１】
以上のアルゴリズムを図３の処理ステップと対応付けて説明すると、次のようになる。まず、処理装置は、集合ＶよりグラフＧを生成し、各頂点に乱数を割り当てる（ステップＳ１３）。次に、２つの頂点の乱数の値を比較して、大きな値に対応する頂点を除き、残った頂点の集合を独立集合Ｉとする（ステップＳ１４）。
【００４２】
次に、集合Ｖ、Ｉ、およびＦを更新する（ステップＳ１５）。ここでは、集合Ｖから、独立集合Ｉおよび直接隣接集合Ｎ（Ｉ）を除いて、新たな集合Ｖを生成する。また、Ｎ（Ｉ）を直接隣接集合の和集合Ｎに加え、独立集合Ｉを最大独立集合Ｆに加える。直接隣接集合Ｎは、グラフのすべての頂点から最大独立集合Ｆの要素を除くことによっても、求めることができる。
【００４３】
次に、集合の更新を終了するか否かを判定する（ステップＳ１６）。ここでは、独立集合ＩがＮＵＬＬになるか、または、ＶまたはＥ（残りの頂点またはエッジ）が空になったときに、更新を終了するものとする。
【００４４】
更新を終了しない場合は、ステップＳ１３以降の処理を繰り返し、更新が終了すると、次に、集合Ｖを最大独立集合Ｆに加えて、最終的な最大独立集合Ｆを決定する（ステップＳ１７）。このとき、最大独立フラグ配列においてフラグがｏｎになっている頂点が、集合Ｆの変数番号を表している。また、直接隣接集合Ｎに属するためにグラフから除かれた頂点であることを示すフラグを格納する１次元配列（直接隣接フラグ配列）も用意され、対応する頂点のフラグがｏｎに設定される。
【００４５】
次に、処理装置は、得られた最大独立集合の変数番号の変換表を作成する（ステップＳ１８）。ここでは、グラフの頂点から最大独立集合に属する変数番号を識別するために、最大独立フラグ配列を参照し、最終的にフラグがｏｎの要素に、最大独立集合に属する頂点の通し番号を振り直す。そして、その番号を、最大独立フラグ配列のフラグの代わりに設定する。これにより、最大独立フラグ配列を利用して、元の変数番号の集合から最大独立集合の新たな変数番号への変換表が作成される。
【００４６】
次に、直接隣接集合の変数番号の変換表を作成する（ステップＳ１９）。この場合、最大独立集合と同様にして、直接隣接フラグ配列を参照し、フラグがｏｎの要素に、直接隣接集合に属する頂点の通し番号を振り直す。これにより、直接隣接フラグ配列を利用して、元の変数番号の集合から直接隣接集合の新たな変数番号への変換表が作成される。
【００４７】
次に、得られた変換表を利用して、ＡFFとＡCFを分離する（図４のステップＳ２０）。このとき、Ａの行ベクトルの要素に関して、変換表の新たな変数番号を参照し、最大独立集合に属する頂点だけを取り出して、それらに対応する行ベクトルの要素を左に詰めて配列ＡＥに格納する。また、これに応じて配列ＩＣＯＬを書き直す。そして、ＡＥとＩＣＯＬの行ベクトルをＡFFとＡCFに振り分けて、ＡFFとＡCFを分離する。
【００４８】
次に、ＡFFの近似逆行列ＦFFを反復法で求める（ステップＳ２１）。ＡFFは対角優位なので、その対角要素からなる対角行列ＤFFが存在する。したがって、ＤFFの逆行列ＤFF -1を初期解としてＭFF -1を反復法により求めれば、この計算は収束する。このとき、｜Ｉ−ＭFF -1ＡFF｜を反復的に計算して期待する近似精度でＭFF -1を求め、得られたＭFF -1をＦFFとする。具体的には、ＭFF,1 -1＝ＤFF -1とおいて、再帰的に以下の計算を行なう。
ＭFF,p -1＝ＭFF,p-1 -1＋（Ｉ−ＭFF,p-1 -1ＡFF）ＤFF -1 （６）
ρp ＝｜Ｉ−ＭFF,p -1ＡFF｜ （７）
そして、τF を所定の定数として、ρp ＜τF を満たしたところまでの近似を採用する。（６）式には、エルパック形式で格納されたスパース行列同士の行列積が現れるが、この計算は以下の手順で行われる。
【００４９】
一般に、Ａ＝Ａ＋αＢ×Ｃの形の行列演算では、まず、メモリの一時領域でｔ＝αＢ×Ｃを計算してから、Ａ＝Ａ＋ｔを計算する。ｔ＝αＢ×Ｃの計算では、エルパック形式で格納されたＢおよびＣの各々の行ベクトルの非ゼロ要素の最大数を、それぞれ、ｎｒｂおよびｎｒｃとして、一時的に第２次元の大きさがｎｒｂ×ｎｒｃの領域を確保する。次に、Ｃの列ベクトルに対してＢの各要素を掛ける。α＝１の場合の手順は、以下の通りである。
（１）ｔの配列ｔｍｐを計算する。
【００５０】
ｉｉ＝１
ｄｏ ｉ＝１，ｎｒｃ
ｄｏ ｋ＝１，ｎｒｂ
ｄｏ ｊ＝１，ｎ
ｔｍｐ（ｊ，ｉｉ）＝ｂ（ｊ，ｋ）＊ｃ（ｂｉｃｏｌ（ｊ，ｋ），ｉ）
ｔｍｐｉｃｏｌ（ｊ，ｉｉ）＝ｃｉｃｏｌ（ｂｉｃｏｌ（ｊ，ｋ），ｉ）
ｅｎｄｄｏ
ｉｉ＝ｉｉ＋１
ｅｎｄｄｏ
ｅｎｄｄｏ
（２）ｔｍｐの行ベクトルの要素を、それらが属する列ベクトルの番号順にソートする。
（３）Ａの行ベクトルも、同様にしてソートする。
（４）Ａ＝Ａ＋ｔｍｐを計算する。ここでは、行ベクトル毎に以下の処理を行う。まず、Ａ（ｉ，ｊ）のＩＣＯＬ（ｉ，ｊ）に一致するｔｍｐｉｃｏｌ（ｉ，ｋ）を、ｋを増加させながら探す。一致する要素が見つかれば、Ａ（ｉ，ｊ）にｔｍｐ（ｉ，ｋ）を加える。また、ｔｍｐｉｃｏｌ（ｉ，ｋ）と一致するＩＣＯＬ（ｉ，ｊ）がなければ、ｔｍｐ（ｉ，ｋ）およびｔｍｐｉｃｏｌ（ｉ，ｋ）を、ＡおよびＩＣＯＬの第ｉ行に新しい要素として追加する。
【００５１】
次に、前処理の計算で用いられるＡCFＦFFを計算し（ステップＳ２２）、ＡFCとＡCCを分離する（ステップＳ２３）。このとき、ステップＳ２０の処理と同様に、Ａの行ベクトルの要素に関して、変換表の新たな変数番号を参照し、直接隣接集合に属する頂点だけを取り出し、配列ＩＣＯＬを書き直す。そして、行ベクトルをＡFCとＡCCに振り分ける。
【００５２】
次に、ＡC ＝ＡCC−ＡCFＦFFＡFCを計算する（ステップＳ２４）。ここでは、ステップＳ２１の行列積の計算と同様にして、ＡCFＦFFＡFCをｔｍｐに格納し、ＡCFＦFFＡFCをＡCCから差し引く。このとき、ＡC の各要素の絶対値を判定しながら、与えられた値より小さな要素を落として、ＭC ＝ＡC −ＲC によりＡC を近似し、ＡC のスパース性を保つ。そして、ＭC をＡC の計算結果として採用する。
【００５３】
このとき、例えば、τを所定の定数として、（ＡC の要素の絶対値）＜τ×（要素の属する行ベクトルの各要素の絶対値の和）という条件を満たすか否かを判定し、この条件を満たす要素を落とすことにする。具体的には、ＲC を決定するための条件は、｜ＭC -1ＲC ｜1 ＜τで与えられる。ただし、ｎ×ｎの行列Ｂ＝（ｂi,j ）に対して、｜Ｂ｜1 は次式で定義される。
【００５４】
【数４】
ここで、未知行列ＭC を以下の対角行列で近似して考えることにする。
【００５６】
【数５】
ここで、ｓｉｇｎ（ｘ）は、ｘ≧０のとき１となり、ｘ＜０のとき−１となる。このような近似によれば、ＡC の要素ａi,j が次式を満たすときに、その要素を落とすことになる。
｜ａi,j ｜＜τ×ｍi,i （１０）
実際には、ＡC の行ベクトルの各要素の絶対値の和を求めて、その行の各要素に関して、左から順に（１０）式の条件を判定し、条件を満たせばその要素を削除する。そして、次の要素を詰めて格納する。
【００５８】
次に、得られたＡC が分解の最も粗いレベル（最終レベル）に対応するか否かを判定する（ステップＳ２５）。最終レベルに対応しなければ、得られたＡC をＡとして用いて、図３のステップＳ１１以降の処理を繰り返す。そして、最終レベルに達したら、行列の分解を終了し、その結果を利用して前処理の計算を行い（ステップＳ２６）、処理を終了する。前処理の計算においては、上述したように、各レベルの不完全ブロック分解を利用して、行列ベクトル積Ｍ-1ｑが求められる。
【００５９】
次に、図６から図９までを参照しながら、上述した前処理を並列に行なう構成について説明する。
図６は、対称型マルチプロセッサ（Symmetric Multi-Processors，ＳＭＰ）システム（並列計算機）の構成図である。図６のマルチプロセッサシステムは、複数のプロセッサ１１、入力装置１２、出力装置１３、および共有メモリ１４を備え、それらはシステムバス１５により接続されている。入力装置１２は、ユーザからの指示や情報の入力に用いられ、出力装置１３は、ユーザへの問い合わせや処理結果の出力に用いられる。
【００６０】
プロセッサ１１は、キャッシュ付きのＣＰＵ（中央処理装置）に対応し、共有メモリ１４を介して互いに結合されている。この場合、並列処理のアプリケーションプログラムは、スレッド（thread）と呼ばれる処理単位を利用して、ＳＭＰの並列拡張言語であるＯｐｅｎＭＰにより記述される。
【００６１】
ＯｐｅｎＭＰでは、変数にＳＨＡＲＥＤおよびＰＲＩＶＡＴＥの属性を付与することができる。このうち、ＳＨＡＲＥＤ属性を持つ変数は、スレッド間で共通に参照・書き込みを行うことができる。また、ＯｐｅｎＭＰでは、ＢＡＲＲＩＥＲ同期を利用して、各スレッドの番号を判定しながら、並列処理を制御することができる。
【００６２】
このシステムでは、処理対象の頂点番号がプロセッサの台数で分割される。そして、図７に示すように、配列ＡＥおよびＩＣＯＬにエルパック形式で格納された行列Ａにおいて、各プロセッサに割り当てられた番号に関する要素の計算を、各スレッドが担当する。ここでは、互いに異なるプロセッサ上に生成されたスレッドＴ１、Ｔ２、Ｔ３、およびＴ４により、行列Ａが分割されている。
【００６３】
ＡFF、ＡCF、ＡFC、およびＡCCの各スレッドの担当部分は、スレッド毎に連続な共用領域を作業領域として用いて格納される。このとき、図８に示すように、行列の非ゼロ要素が存在する部分のバンド幅（上バンド幅および下バンド幅）が求められ、各スレッドの担当部分のバンド幅が決められる。
【００６４】
また、図３のステップＳ１３で各頂点に割り当てられる乱数を格納する配列、最大独立フラグ配列、直接隣接フラグ配列、および変数番号の変換表等の格納領域は、図９に示すように設定される。図９では、各スレッドが担当する部分の格納領域２２の左側に、下バンド幅分の余分な領域（オーバラップ領域）２１が設けられ、右側に上バンド幅分のオーバラップ領域２３が設けられている。
【００６５】
ステップＳ１３〜Ｓ１７の処理においては、最大独立集合を決定するアルゴリズムを以下の方法で並列化する。ここでは、フラグの“ｏｎ”と“ｏｆｆ”をそれぞれ整数の“１”と“０”で表し、組み込み関数ＭＩＮを利用して設定する。最大独立フラグ配列の初期値としては、直接隣接集合の頂点以外をｏｎとし、直接隣接集合として除かれた頂点をｏｆｆとする。
【００６６】
まず、各スレッドは、乱数を生成して、グラフを構成する頂点のオーバラップ領域以外の領域（各スレッドが担当する領域）に格納する。その後、オーバラップ領域に、その部分を担当する他のスレッドの領域から乱数をコピーして、値を設定し直す。これにより、各スレッドで用いる頂点の乱数が、連続領域でローカルに保持される。
【００６７】
次に、各スレッドは、担当するグラフを表現するエッジの頂点対をすべて処理する。処理対象の頂点は、各スレッドが担当する部分の頂点番号に上バンド幅分および下バンド幅分の頂点番号を加えた範囲に存在する。そして、頂点対に対応する２つの乱数を比較し、最大独立フラグ配列において、値が大きい方の頂点のフラグをｏｆｆにする。
【００６８】
各スレッドでの設定が終了すると、各スレッドは、オーバラップ領域の頂点を担当する他のスレッドからフラグ情報をコピーして、最大独立フラグ配列のオーバラップ領域に格納する。このとき、関数ＭＩＮを用いてスレッド間で整合性をとり、他のスレッドでｏｆｆになっている部分はｏｆｆにする。
【００６９】
次に、各スレッドは、もう一度グラフのエッジを走査し、一方の頂点が最大独立集合に属する場合、直接隣接フラグ配列において、もう一方の頂点のフラグをｏｎにする。また、最大独立フラグ配列と同様にして、スレッド間で整合性をとるように、直接隣接フラグ配列のオーバラップ領域の値を設定する。
【００７０】
次に、各スレッドは、グラフを再度走査し、最大独立集合と直接隣接集合に属さないエッジのみを、グラフを構成するエッジとして残す。このとき、グラフを構成するエッジだけを詰め、最大独立集合の要素数をカウントして、それがゼロか否かをチェックする。また、頂点およびエッジの数をカウントして、それらがゼロか否かをチェックする。最大独立集合の要素数、頂点およびエッジの数が共にゼロでなければ、新たなグラフに関して、乱数の生成以降の処理を繰り返す。
【００７１】
そして、最大独立集合の要素数またはエッジの数がゼロになると、各スレッドは、繰り返しを終了し、直接隣接フラグ配列において、エッジを構成する頂点のフラグをｏｎにする。そして、スレッド間で格納領域の整合性をとる。
【００７２】
次に、ステップＳ１８において、アプリケーションは、各スレッドが担当する頂点（オーバラップ領域は除く）のうち、最大独立集合に属するものの数をカウントする。そして、スレッドの総数の大きさの１次元配列（共有領域）に各スレッドのカウント値を格納し、この情報を用いて、各スレッドの最大独立集合の頂点に、全体での通し番号を振る。このとき、通し番号は、正整数で１から順に振られる。
【００７３】
次に、各スレッドは、最大独立フラグ配列のｏｎの代わりに、その頂点に対応する通し番号を書き込み、変数番号の変換表を作成する。この変換表において、正整数に対応する変数番号は、最大独立集合に属する頂点に対応し、０に対応する変数番号は、直接隣接集合に属する頂点に対応する。そして、オーバラップ領域に対応する他のスレッドから、変換表の一部をコピーする。ステップＳ１９においても、同様にして、直接隣接集合の変数番号の変換表が作成される。
【００７４】
次に、図４のステップＳ２０において、各スレッドは、並列に作業領域を確保し、最大独立集合の変換表を用いてＡFFとＡCFを分離する。この場合、Ａの要素ａi,j のうち、ｊが最大独立集合に属するもののみを、ＡＥ内で左に詰めて格納し、ＩＣＯＬの格納値をその頂点の新たな変数番号に書き換える。
【００７５】
次に、各行ベクトルの非ゼロ要素の個数の最大値Ｎｒを求め、各スレッドが担当する最大独立集合および直接隣接集合の要素数とＮｒに基づいて、各スレッドの作業領域にＡFFとＡCFの格納領域（ＡＥおよびＩＣＯＬ）を確保する。そして、ＡのＡＥおよびＩＣＯＬの行ベクトルをＡFFとＡCFの格納領域に振り分ける。ステップＳ２３においても、同様にして、直接隣接集合の変換表を用いてＡFCとＡCCが分離される。
【００７６】
また、ステップＳ２１およびＳ２４においては、スパース行列積の計算が並列に行われる。各スレッドは、図７および図８と同様に、処理対象の各行列の対応部分を保持している。
【００７７】
上述したｔｍｐの計算において、最も内側のｄｏループの処理は、各スレッドに割り当てられた部分に対してのみ行われる。また、必要なＣの列ベクトルに関しても、上バンド幅と下バンド幅のオーバラップ領域を有する作業領域を用意し、各スレッドが保持する部分と他のスレッドが保持する部分をその作業領域にコピーしてから、各スレッドで並列に計算する。
【００７８】
また、ステップＳ２６のｐ＝Ｍ-1ｑの並列計算において、アプリケーションは、ベクトルｑF およびｑC の各々に対して、各スレッドが担当する部分の格納領域と上下バンド幅のオーバラップ領域を用意し、変換表を利用して、ベクトルｑからｑF およびｑC を分離する。そして、各スレッドは、割り当てられた行列部分を用いて、行列ベクトル積を並列に計算する。このとき、他のスレッドの処理結果を参照するために、オーバラップ領域を用いる。
【００７９】
ところで、本実施形態の処理装置は、並列計算機に限らず、図１０に示すような任意の情報処理装置（コンピュータ）を用いて構成することもできる。図１０の情報処理装置は、ＣＰＵ（中央処理装置）３１、メモリ３２、入力装置３３、出力装置３４、外部記憶装置３５、媒体駆動装置３６、およびネットワーク接続装置３７を備え、それらはバス３８により互いに接続されている。
【００８０】
メモリ３２は、例えば、ＲＯＭ（read only memory）、ＲＡＭ（random access memory）等を含み、処理に用いられるプログラムとデータを格納する。ＣＰＵ３１は、メモリ３２を利用してプログラムを実行することにより、必要な処理を行う。上述した連立１次方程式の反復法および前処理のプログラムは、インストラクションの集合として、メモリ３２の特定のプログラムコードセグメントに格納される。
【００８１】
入力装置３３は、例えば、キーボード、ポインティングデバイス、タッチパネル等であり、ユーザからの指示や情報の入力に用いられる。出力装置３４は、例えば、ディスプレイ、プリンタ、スピーカ等であり、ユーザへの問い合わせや処理結果の出力に用いられる。
【００８２】
外部記憶装置３５は、例えば、磁気ディスク装置、光ディスク装置、光磁気ディスク（magneto-optical disk）装置、テープ装置等である。情報処理装置は、この外部記憶装置３５に、上述のプログラムとデータを保存しておき、必要に応じて、それらをメモリ３２にロードして使用する。
【００８３】
媒体駆動装置３６は、可搬記録媒体３９を駆動し、その記録内容にアクセスする。可搬記録媒体３９としては、メモリカード、フロッピーディスク、ＣＤ−ＲＯＭ（compact disk read only memory ）、光ディスク、光磁気ディスク等、任意のコンピュータ読み取り可能な記録媒体が用いられる。ユーザは、この可搬記録媒体３９に上述のプログラムとデータを格納しておき、必要に応じて、それらをメモリ３２にロードして使用する。
【００８４】
ネットワーク接続装置３７は、ＬＡＮ（local area network）等の任意の通信ネットワークに接続され、通信に伴うデータ変換を行う。また、情報処理装置は、上述のプログラムとデータをネットワーク接続装置３７を介して、サーバ等の他の装置から受け取り、必要に応じて、それらをメモリ３２にロードして使用する。
【００８５】
図１１は、図１０の情報処理装置にプログラムとデータを供給することのできるコンピュータ読み取り可能な記録媒体を示している。可搬記録媒体３９やサーバ４０のデータベース４１に保存されたプログラムとデータは、メモリ３２にロードされる。そして、ＣＰＵ３１は、そのデータを用いてそのプログラムを実行し、必要な処理を行う。このとき、サーバ４０は、プログラムとデータを伝送する伝搬信号を生成し、ネットワーク上の任意の伝送媒体を介して、情報処理装置に送信する。
（付記１） 連立一次方程式を解く反復法のための前処理を、マルチレベル不完全ブロック分解により行なう処理装置であって、
前記不完全ブロック分解のあるレベルにおける係数行列を格納する係数格納手段と、
前記係数格納手段に格納された係数行列の要素のうち、取り除かれる変数の変数番号を行番号および列番号とする要素からなるブロック行列が対角優位になるように、該取り除かれる変数の変数番号の集合を決定する決定手段と、
決定された変数番号の集合に関する情報を格納する集合格納手段と、
前記集合格納手段に格納された情報に基づいて、前記ブロック行列の近似逆行列を求める逆行列手段と、
前記係数格納手段に格納された係数行列の要素のうち、残される変数の変数番号を行番号および列番号とする要素からなるブロック行列から、粗いレベルの係数行列を求める係数行列手段と、
次の不完全ブロック分解を行うために、前記粗いレベルの係数行列を前記係数格納手段に格納する分解手段と、
前記近似逆行列を用いて、前記反復法における行列ベクトル積を計算する計算手段と
を備えることを特徴とする処理装置。
（付記２） 前記決定手段は、前記係数格納手段に格納された係数行列の要素の添え字の集合の部分集合である直接隣接集合の初期値を設定し、該直接隣接集合の初期値と乱数を用いて添え字の最大独立集合を求めることで、前記変数番号の集合を決定することを特徴とする付記１記載の処理装置。
（付記３） 前記逆行列手段は、与えられた近似精度で反復計算により前記近似逆行列を求めることを特徴とする付記１記載の処理装置。
（付記４） 前記係数行列手段は、前記粗いレベルの係数行列の非ゼロ要素を削減して、該粗いレベルの係数行列のスパース性を保持することを特徴とする付記１記載の処理装置。
（付記５） 複数のプロセッサを有し、連立一次方程式を解く反復法のための前処理を、マルチレベル不完全ブロック分解により行なう並列計算機であって、
前記複数のプロセッサが前記不完全ブロック分解のあるレベルにおける係数行列を並列に処理できるように、該係数行列を分割して格納する係数格納手段と、前記係数格納手段に格納された係数行列の要素のうち、取り除かれる変数の変数番号を行番号および列番号とする要素からなるブロック行列が対角優位になるように、該取り除かれる変数の変数番号の集合を並列に決定する決定手段と、
前記複数のプロセッサが決定された変数番号の集合に関する情報を並列に処理できるように、該変数番号の集合に関する情報を分割して格納する集合格納手段と、
前記集合格納手段に格納された情報に基づいて、前記ブロック行列の近似逆行列を並列に求める逆行列手段と、
前記係数格納手段に格納された係数行列の要素のうち、残される変数の変数番号を行番号および列番号とする要素からなるブロック行列から、粗いレベルの係数行列を並列に求める係数行列手段と、
次の不完全ブロック分解を行うために、前記粗いレベルの係数行列を前記係数格納手段に分割して格納する分解手段と、
前記近似逆行列を用いて、前記反復法における行列ベクトル積を並列に計算する計算手段と
を備えることを特徴とする並列計算機。
（付記６） 連立一次方程式を解く反復法のための前処理を、マルチレベル不完全ブロック分解により行なうコンピュータのためのプログラムを記録した記録媒体であって、該プログラムは、
前記不完全ブロック分解のあるレベルにおける係数行列を格納し、
格納された係数行列の要素のうち、取り除かれる変数の変数番号を行番号および列番号とする要素からなるブロック行列が対角優位になるように、該取り除かれる変数の変数番号の集合を決定し、
決定された変数番号の集合に関する情報を格納し、
格納された情報に基づいて、前記ブロック行列の近似逆行列を求め、
前記格納された係数行列の要素のうち、残される変数の変数番号を行番号および列番号とする要素からなるブロック行列から、粗いレベルの係数行列を求め、次の不完全ブロック分解を行うために、前記粗いレベルの係数行列を前記あるレベルにおける係数行列として格納し、
前記近似逆行列を用いて、前記反復法における行列ベクトル積を計算する
処理を前記コンピュータに実行させることを特徴とするコンピュータ読み取り可能な記録媒体。
（付記７） 連立一次方程式を解く反復法のための前処理を、マルチレベル不完全ブロック分解により行なうコンピュータにプログラムを伝送する伝搬信号であって、該プログラムは、
前記不完全ブロック分解のあるレベルにおける係数行列を格納し、
格納された係数行列の要素のうち、取り除かれる変数の変数番号を行番号および列番号とする要素からなるブロック行列が対角優位になるように、該取り除かれる変数の変数番号の集合を決定し、
決定された変数番号の集合に関する情報を格納し、
格納された情報に基づいて、前記ブロック行列の近似逆行列を求め、
前記格納された係数行列の要素のうち、残される変数の変数番号を行番号および列番号とする要素からなるブロック行列から、粗いレベルの係数行列を求め、
次の不完全ブロック分解を行うために、前記粗いレベルの係数行列を前記あるレベルにおける係数行列として格納し、
前記近似逆行列を用いて、前記反復法における行列ベクトル積を計算する
処理を前記コンピュータに実行させることを特徴とする伝搬信号。
【００８６】
【発明の効果】
本発明によれば、係数行列がＭ行列でないような連立一次方程式の反復法において、計算の収束を速めるような前処理を行うことが可能になる。したがって、従来の前処理では収束しなかったような問題を、効率良く解くことができる。
【図面の簡単な説明】
【図１】本発明の処理装置の原理図である。
【図２】マルチレベル不完全ブロック分解を示す図である。
【図３】前処理のフローチャート（その１）である。
【図４】前処理のフローチャート（その２）である。
【図５】エルパック形式格納法を示す図である。
【図６】マルチプロセッサシステムの構成図である。
【図７】配列の分割を示す図である。
【図８】バンド幅を示す図である。
【図９】格納領域を示す図である。
【図１０】情報処理装置の構成図である。
【図１１】記録媒体を示す図である。
【符号の説明】
１ 係数格納手段
２ 決定手段
３ 集合格納手段
４ 逆行列手段
５ 係数行列手段
６ 分解手段
７ 計算手段
１１ プロセッサ
１２、３３ 入力装置
１３、３４ 出力装置
１４ 共有メモリ
１５ システムバス
２１、２３ オーバラップ領域
２２ 格納領域
３１ ＣＰＵ
３２ メモリ
３５ 外部記憶装置
３６ 媒体駆動装置
３７ ネットワーク接続装置
３８ バス
３９ 可搬記録媒体
４０ サーバ
４１ データベース
Claims (4)
- 係数行列がＭ行列でない連立一次方程式を解く反復法のための前処理を、マルチレベル不完全ブロック分解により行なう処理装置であって、
前記不完全ブロック分解のあるレベルにおける係数行列を格納する係数格納手段と、
前記係数格納手段に格納された係数行列の要素のうち、取り除かれる変数の変数番号を行番号および列番号とする要素からなるブロック行列が対角優位になるように、該取り除かれる変数の変数番号の集合を決定する決定手段と、
決定された変数番号の集合に関する情報を格納する集合格納手段と、
前記集合格納手段に格納された情報に基づいて、前記ブロック行列の近似逆行列を求める逆行列手段と、
前記係数格納手段に格納された係数行列の要素のうち、残される変数の変数番号を行番号および列番号とする要素からなるブロック行列から、次のレベルの係数行列を求める係数行列手段と、
次の不完全ブロック分解を行うために、前記次のレベルの係数行列を前記係数格納手段に格納する分解手段と、
前記近似逆行列を用いて、前記反復法における行列ベクトル積を計算する計算手段とを備え、
前記決定手段は、前記係数格納手段に格納された係数行列の要素のうち、対角優位な各要素の添え字に乱数を割り当て、２つの添え字に割り当てられた乱数を比較して大きい方に対応する添え字を取り除く処理を繰り返すことで、残された添え字からなる最大独立集合を求め、該最大独立集合を前記変数番号の集合とすることを特徴とする処理装置。 - 前記係数行列手段は、前記次のレベルの係数行列の非ゼロ要素を削減して、該次のレベルの係数行列のスパース性を保持することを特徴とする請求項１記載の処理装置。
- 複数のプロセッサを有し、係数行列がＭ行列でない連立一次方程式を解く反復法のための前処理を、マルチレベル不完全ブロック分解により行なう並列計算機であって、
前記複数のプロセッサが前記不完全ブロック分解のあるレベルにおける係数行列を並列に処理できるように、該係数行列を分割して格納する係数格納手段と、
前記係数格納手段に格納された係数行列の要素のうち、取り除かれる変数の変数番号を行番号および列番号とする要素からなるブロック行列が対角優位になるように、該取り除かれる変数の変数番号の集合を並列に決定する決定手段と、
前記複数のプロセッサが決定された変数番号の集合に関する情報を並列に処理できるように、該変数番号の集合に関する情報を分割して格納する集合格納手段と、
前記集合格納手段に格納された情報に基づいて、前記ブロック行列の近似逆行列を並列に求める逆行列手段と、
前記係数格納手段に格納された係数行列の要素のうち、残される変数の変数番号を行番号および列番号とする要素からなるブロック行列から、次のレベルの係数行列を並列に求める係数行列手段と、
次の不完全ブロック分解を行うために、前記次のレベルの係数行列を前記係数格納手段に分割して格納する分解手段と、
前記近似逆行列を用いて、前記反復法における行列ベクトル積を並列に計算する計算手段とを備え、
前記決定手段は、前記係数格納手段に格納された係数行列の要素のうち、対角優位な各要素の添え字に乱数を割り当て、２つの添え字に割り当てられた乱数を比較して大きい方に対応する添え字を取り除く処理を繰り返すことで、残された添え字からなる最大独立集合を求め、該最大独立集合を前記変数番号の集合とすることを特徴とする並列計算機。 - 係数行列がＭ行列でない連立一次方程式を解く反復法のための前処理を、マルチレベル不完全ブロック分解により行なうコンピュータのためのプログラムを記録した記録媒体であって、該プログラムは、
前記不完全ブロック分解のあるレベルにおける係数行列を係数格納手段に格納し、
前記係数格納手段に格納された係数行列の要素のうち、取り除かれる変数の変数番号を行番号および列番号とする要素からなるブロック行列が対角優位になるように、該取り除かれる変数の変数番号の集合を決定し、
決定された変数番号の集合に関する情報を集合格納手段に格納し、
前記集合格納手段に格納された情報に基づいて、前記ブロック行列の近似逆行列を求め、
前記係数格納手段に格納された係数行列の要素のうち、残される変数の変数番号を行番号および列番号とする要素からなるブロック行列から、次のレベルの係数行列を求め、
次の不完全ブロック分解を行うために、前記次のレベルの係数行列を前記あるレベルにおける係数行列として前記係数格納手段に格納し、
前記近似逆行列を用いて、前記反復法における行列ベクトル積を計算する
処理を前記コンピュータに実行させ、
前記コンピュータは、前記変数番号の集合を決定する処理において、前記係数格納手段に格納された係数行列の要素のうち、対角優位な各要素の添え字に乱数を割り当て、２つの添え字に割り当てられた乱数を比較して大きい方に対応する添え字を取り除く処理を繰り返すことで、残された添え字からなる最大独立集合を求め、該最大独立集合を該変数番号の集合とすることを特徴とするコンピュータ読み取り可能な記録媒体。
Priority Applications (2)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| JP2000354179A JP3809062B2 (ja) | 2000-11-21 | 2000-11-21 | マルチレベル不完全ブロック分解による前処理を行う処理装置 |
| US09/891,203 US6799194B2 (en) | 2000-11-21 | 2001-06-26 | Processing apparatus for performing preconditioning process through multilevel block incomplete factorization |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| JP2000354179A JP3809062B2 (ja) | 2000-11-21 | 2000-11-21 | マルチレベル不完全ブロック分解による前処理を行う処理装置 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| JP2002157237A JP2002157237A (ja) | 2002-05-31 |
| JP3809062B2 true JP3809062B2 (ja) | 2006-08-16 |
Family
ID=18826823
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| JP2000354179A Expired - Fee Related JP3809062B2 (ja) | 2000-11-21 | 2000-11-21 | マルチレベル不完全ブロック分解による前処理を行う処理装置 |
Country Status (2)
| Country | Link |
|---|---|
| US (1) | US6799194B2 (ja) |
| JP (1) | JP3809062B2 (ja) |
Families Citing this family (8)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US7069556B2 (en) * | 2001-09-27 | 2006-06-27 | Intel Corporation | Method and apparatus for implementing a parallel construct comprised of a single task |
| US7493603B2 (en) * | 2002-10-15 | 2009-02-17 | International Business Machines Corporation | Annotated automaton encoding of XML schema for high performance schema validation |
| JP4167889B2 (ja) * | 2002-12-06 | 2008-10-22 | 株式会社日立製作所 | データ変換方法、および、そのための計算機システム |
| US7437374B2 (en) * | 2004-02-10 | 2008-10-14 | International Business Machines Corporation | Efficient XML schema validation of XML fragments using annotated automaton encoding |
| US8249451B2 (en) * | 2007-08-16 | 2012-08-21 | Futurewei Technologies, Inc. | Methods for characterizing optical switches and multiplexers/demultiplexers |
| ATE518171T1 (de) * | 2008-05-29 | 2011-08-15 | Dspace Gmbh | Verfahren und simulator zur echtzeitberechnung der zustandsgrössen eines prozessmodells |
| EP2350915A4 (en) * | 2008-09-30 | 2013-06-05 | Exxonmobil Upstream Res Co | METHOD FOR SOLVING STORAGE SIMULATION MATRIX COMPENSATION USING PARALLEL INCOMPLETE MULTILEVEL FACTORIZATION |
| CA2933809C (en) | 2014-01-31 | 2017-07-18 | Landmark Graphics Corporation | Flexible block ilu factorization |
Family Cites Families (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP2956800B2 (ja) * | 1991-09-19 | 1999-10-04 | 株式会社日立製作所 | 連立一次方程式に関する計算装置 |
| JP2959525B2 (ja) * | 1997-06-02 | 1999-10-06 | 日本電気株式会社 | データ処理装置および方法、情報記憶媒体 |
| TW388921B (en) * | 1997-11-28 | 2000-05-01 | Nippon Electric Co | Semiconductor process device simulation method and storage medium storing simulation program |
| JP3639206B2 (ja) * | 2000-11-24 | 2005-04-20 | 富士通株式会社 | 共有メモリ型スカラ並列計算機における並列行列処理方法、及び記録媒体 |
-
2000
- 2000-11-21 JP JP2000354179A patent/JP3809062B2/ja not_active Expired - Fee Related
-
2001
- 2001-06-26 US US09/891,203 patent/US6799194B2/en not_active Expired - Fee Related
Also Published As
| Publication number | Publication date |
|---|---|
| US6799194B2 (en) | 2004-09-28 |
| JP2002157237A (ja) | 2002-05-31 |
| US20020099748A1 (en) | 2002-07-25 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Altschuler et al. | Near-linear time approximation algorithms for optimal transport via Sinkhorn iteration | |
| Pednault et al. | Breaking the 49-qubit barrier in the simulation of quantum circuits | |
| US10013477B2 (en) | Accelerated discrete distribution clustering under wasserstein distance | |
| Fu et al. | Simple hardware-efficient long convolutions for sequence modeling | |
| Krishnan et al. | Efficient preconditioning of laplacian matrices for computer graphics | |
| Kannan et al. | MPI-FAUN: An MPI-based framework for alternating-updating nonnegative matrix factorization | |
| JP2956800B2 (ja) | 連立一次方程式に関する計算装置 | |
| Ravi et al. | Large scale distributed semi-supervised learning using streaming approximation | |
| US11971906B2 (en) | Clustering apparatus, clustering method, program and data structure | |
| JP3809062B2 (ja) | マルチレベル不完全ブロック分解による前処理を行う処理装置 | |
| Takezawa et al. | Supervised tree-wasserstein distance | |
| Lahn et al. | A combinatorial algorithm for approximating the optimal transport in the parallel and mpc settings | |
| Singh et al. | Exploiting sparsity in pruned neural networks to optimize large model training | |
| Li et al. | An alternating nonmonotone projected Barzilai–Borwein algorithm of nonnegative factorization of big matrices | |
| US11042604B2 (en) | Assessing distances between pairs of histograms based on relaxed flow constraints | |
| Zeng et al. | Slice-based online convolutional dictionary learning | |
| Trifunovic et al. | A parallel algorithm for multilevel k-way hypergraph partitioning | |
| Dong et al. | An optimization method for pruning rates of each layer in CNN based on the GA-SMSM | |
| Wang et al. | Learning Large-Scale MTP $ _2 $ Gaussian Graphical Models via Bridge-Block Decomposition | |
| Felip et al. | Tree pyramidal adaptive importance sampling | |
| Meng et al. | Parallel edge-based visual assessment of cluster tendency on GPU | |
| Duchi et al. | Conic descent and its application to memory-efficient optimization over positive semidefinite matrices | |
| Zdunek et al. | Image completion with approximate convex hull tensor decomposition | |
| Li et al. | CUSNTF: A scalable sparse non-negative tensor factorization model for large-scale industrial applications on multi-GPU | |
| Shivdikar et al. | Speeding up dnns using hpl based fine-grained tiling for distributed multi-gpu training |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| A977 | Report on retrieval |
Free format text: JAPANESE INTERMEDIATE CODE: A971007 Effective date: 20051101 |
|
| A131 | Notification of reasons for refusal |
Free format text: JAPANESE INTERMEDIATE CODE: A131 Effective date: 20051108 |
|
| A521 | Request for written amendment filed |
Free format text: JAPANESE INTERMEDIATE CODE: A523 Effective date: 20051227 |
|
| TRDD | Decision of grant or rejection written | ||
| A01 | Written decision to grant a patent or to grant a registration (utility model) |
Free format text: JAPANESE INTERMEDIATE CODE: A01 Effective date: 20060516 |
|
| A61 | First payment of annual fees (during grant procedure) |
Free format text: JAPANESE INTERMEDIATE CODE: A61 Effective date: 20060519 |
|
| R150 | Certificate of patent or registration of utility model |
Free format text: JAPANESE INTERMEDIATE CODE: R150 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20090526 Year of fee payment: 3 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20100526 Year of fee payment: 4 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20100526 Year of fee payment: 4 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20110526 Year of fee payment: 5 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20120526 Year of fee payment: 6 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20130526 Year of fee payment: 7 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20140526 Year of fee payment: 8 |
|
| LAPS | Cancellation because of no payment of annual fees |