(57)【要約】【課題】ガウス・ザイデル法及びＳＯＲ法の少なくともいずれかを並列化して高速に実行することができる並列処理装置等の提供。【解決手段】ガウス・ザイデル法及び逐次加速緩和法（ＳＯＲ法）の少なくともいずれかで行列データを用いて行列連立一次方程式を解く際に、前記行列データの所定行データずつ演算して各行毎に解データを出力する並列処理装置において、前記所定行データの行列データの対角要素の右側を各行毎に並列演算して右側演算結果データを算出する右側演算結果データ算出部と、前記所定行データの前記対角要素の左側と算出済み解データを用いて各行毎に並列演算し対角要素の左側の左側演算結果データを算出する左側演算結果データ算出部と、前記右側演算結果データと前記左側演算結果データとを用いて、所定行データの対角要素の解データを算出する対角要素データ算出部と、を有する並列処理装置である。【選択図】図１０
2018206078.tif 000002
【特許請求の範囲】
【請求項１】
ガウス・ザイデル法及び逐次加速緩和法（ＳＯＲ法）の少なくともいずれかで行列データを用いて行列連立一次方程式を解く際に、前記行列データの所定行データずつ演算して各行毎に解データを出力する並列処理装置において、
前記所定行データの行列データの対角要素の右側を各行毎に並列演算して右側演算結果データを算出する右側演算結果データ算出部と、
前記所定行データの行列データの対角要素の左側と算出済み解データを用い各行毎に並列演算して左側演算結果データを算出する左側演算結果データ算出部と、
前記右側演算結果データと前記左側演算結果データとを用いて、前記所定行データの行列データの対角要素の解データを算出する対角要素解データ算出部と、
を有することを特徴とする並列処理装置。
【請求項２】
前記対角要素の解データを算出する際に、次の所定行データの行列データの対角要素の右側を各行毎に並列演算して右側演算結果データを算出する請求項１に記載の並列処理装置。
【請求項３】
ガウス・ザイデル法及び逐次加速緩和法（ＳＯＲ法）の少なくともいずれかで行列データを用いて行列連立一次方程式を解く際に、前記行列データの所定行データずつ演算して各行毎に解データを出力する並列演算方法において、
前記所定行データの行列データの対角要素の右側を各行毎に並列演算して右側演算結果データを算出し、
前記所定行データの行列データの対角要素の左側と算出済み解データを用い各行毎に並列演算して左側演算結果データを算出し、
前記右側演算結果データと前記左側演算結果データとを用いて、前記所定行データの行列データの対角要素の解データを算出することを特徴とする並列演算方法。
【請求項４】
ガウス・ザイデル法及び逐次加速緩和法（ＳＯＲ法）の少なくともいずれかで行列データを用いて行列連立一次方程式を解く際に、前記行列データの所定行データずつ演算して各行毎に解データを出力する並列演算プログラムにおいて、
前記所定行データの行列データの対角要素の右側を各行毎に並列演算して右側演算結果データを算出し、
前記所定行データの行列データの対角要素の左側と算出済み解データを用い各行毎に並列演算して左側演算結果データを算出し、
前記右側演算結果データと前記左側演算結果データとを用いて、前記所定行データの行列データの対角要素の解データを算出する
処理をコンピュータに実行させることを特徴とする並列演算プログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、並列処理装置、並列演算方法、及び並列演算プログラムに関する。
【背景技術】
【０００２】
流体解析や構造解析などの数値シミュレーションにおいては、疎行列連立一次方程式Ａｘ＝ｂ（Ａ：疎行列、ｘ:解ベクトル、ｂ：部分ベクトル）を解く部分が実行時間の大半を占める場合が多い。そのため、疎行列連立一次方程式を解く部分の高速化は非常に重要である。疎行列とは、行列要素がほとんど零である行列を意味する。
【０００３】
大規模な疎行列連立一次方程式の解法としては、反復法というアルゴリズムが用いられている。反復法では、ｘに適当な初期値を与え、真の解に対して予め与えた許容誤差以下となるまでｘの更新を繰り返すことにより解を求める。解が求まるまでのｘの更新回数を反復回数と呼ぶ。
反復法では、できるだけ反復回数を少なくするために前処理が行われている。有効な前処理の一つとしてマルチグリッド法がある。マルチグリッド法を用いて前処理を行う場合、その中で大きな実行時間を占めるのは、ガウス・ザイデル法やＳＯＲ法（Ｓｕｃｃｅｓｓｉｖｅ  Ｏｖｅｒ−Ｒｅｌａｘａｔｉｏｎ；逐次加速緩和法）による行列計算であることが多い。そのため、効率的に前処理を行うためには、ガウス・ザイデル法やＳＯＲ法の高速化が重要となる。
【０００４】
ここで、図１にガウス・ザイデル法の概要を示す。図１中の（ａ）〜（ｅ）の順に実行が進むにつれて、ベクトルｘの要素をｘ［０］からｘ［ｎ−１］の順に更新する。ｘ［Ａｃｏｌ［ｉ］［ｊ］］（図１中の（ｄ））と、ｘ［ｉ］（図１中の（ｅ））とは依存関係がある。そのため、ｉ＝０行目の計算を行い、次に、ｉ＝１行目の計算を行い、ｉ＝２行目の計算を行い、・・・というように逐次的に実行する。
【０００５】
行列Ａが疎行列である場合、行ｉ１と行ｉ２とは、互いに同じ列位置に非零要素を持たない。即ち、ｘ［Ａｃｏｌ［ｉ１］［ｊ］］とｘ［ｉ２］との間に依存関係が存在せず、かつｘ［Ａｃｏｌ［ｉ２］［ｊ］］とｘ［ｉ１］との間に依存関係が存在しない場合がある。その場合、２つの行ｉ１と行ｉ２の計算を並列に行うことができる。即ち、ｘ［ｉ１］とｘ［ｉ２］の値を同時に更新することができる。
以上のことは、元々のベクトルｘの更新順序、即ち、ｘ［０］からｘ［ｎ−１］まで順に更新することを変えることになり、反復回数が増加してしまう。
【０００６】
一方、ガウス・ザイデル法やＳＯＲ法とは別の連立一次方程式の解法として正方行列を下三角形行列Ｌと上三角形行列Ｕの積に分解するＬＵ分解法がある。このＬＵ分解法で解く際に、Ｋ回目の反復処理を逐次処理Ｐと並列化可能な処理Ｑとに分け、処理ＱをＱ１〜Ｑｎに分解し、Ｑ１とＰ（Ｋ＋Ｍ）を１つのプロセッサで順次実行する。これらの処理と並行して他の処理Ｑ２〜Ｑｎを他のプロセッサで並列実行し、処理順序を保証してプロセッサの待ち状態の時間を軽減することが提案されている（例えば、特許文献１参照）。
【先行技術文献】
【特許文献】
【０００７】
【特許文献１】特開平８−２２７４０５号公報
【発明の概要】
【発明が解決しようとする課題】
【０００８】
しかしながら、前述のＬＵ分解法は、非零要素が多く含まれる密行列を対象とする手法であり、零要素が多く含まれる疎行列を対象とするガウス・ザイデル法やＳＯＲ法とは対象とする行列が異なる。また、前述のＬＵ分解法は直接法であり、反復法であるガウス・ザイデル法やＳＯＲ法とは解を求める手法が異なる。したがって、前述のＬＵ分解法では、演算が進むにつれ演算の並列度が下がり、疎行列に対するガウス・ザイデル法及びＳＯＲ法の少なくともいずれかを高速に実行することはできない。
【課題を解決するための手段】
【０００９】
一つの実施態様では、ガウス・ザイデル法及び逐次加速緩和法（ＳＯＲ法）の少なくともいずれかで行列データを用いて行列連立一次方程式を解く際に、行列データの所定行データずつ演算して各行毎に解データを出力する並列処理装置において、
所定行データの行列データの対角要素の右側を各行毎に並列演算して右側演算結果データを算出する右側演算結果データ算出部と、
所定行データの行列データの対角要素の左側と算出済み解データを用い各行毎に並列演算して左側演算結果データを算出する左側演算結果データ算出部と、
右側演算結果データと左側演算結果データとを用いて、所定行データの行列データの対角要素の解データを算出する対角要素解データ算出部と、を有する並列処理装置である。
【発明の効果】
【００１０】
ガウス・ザイデル法及びＳＯＲ法の少なくともいずれかを並列化して高速に実行することができる並列処理装置、並列演算方法、及び並列演算プログラムを提供することができる。
【図面の簡単な説明】
【００１１】
【図１】図１は、ガウス・ザイデル法の概要の一例を示す図である。
【図２】図２は、ガウス・ザイデル法の概要の他の一例を示す図である。
【図３】図３は、ガウス・ザイデル法による３×３行列の計算の一例を示す図である。
【図４】図４は、従来の並列処理装置全体の処理の一例を示すブロック図である。
【図５】図５は、従来の並列処理装置全体の流れの一例を示すフローチャートである。
【図６】図６は、従来のガウス・ザイデル法の一例を説明するための説明図である。
【図７】図７は、従来のガウス・ザイデル法による処理の流れの一例を示すフローチャートである。
【図８】図８は、従来のガウス・ザイデル法による１６×１６行列の計算の一例を示す図である。
【図９】図９は、本発明の並列処理装置のハードウェア構成の一例を示す図である。
【図１０】図１０は、本発明の並列処理装置全体の処理の一例を示すブロック図である。
【図１１】図１１は、本発明の並列処理装置全体の処理の流れの一例を示すフローチャートである。
【図１２】図１２は、本発明のガウス・ザイデル法の一例を説明するための説明図である。
【図１３】図１３は、本発明のガウス・ザイデル法による処理の流れの一例を示すフローチャートである。
【図１４】図１４は、本発明のガウス・ザイデル法による１６×１６行列の計算の一例を示す図である。
【図１５】図１５は、実施例における疎行列の分割方法の一例を示す図である。
【図１６】図１６は、実施例１における分割した疎行列の実行順序の一例を示す図である。
【図１７】図１７は、実施例２における分割した疎行列の実行順序の一例を示す図である。
【発明を実施するための形態】
【００１２】
以下、本発明の一実施形態について説明するが、本発明は、これらの実施形態に何ら限定されるものではない。
【００１３】
本発明の並列処理装置は、ガウス・ザイデル法及び逐次加速緩和法（ＳＯＲ法）の少なくともいずれかで行列データを用いて行列連立一次方程式を解く際に、行列データの所定行データずつ演算して各行毎に解データを出力する並列処理装置において、右側演算結果データ算出部と、左側演算結果データ算出部と、対角要素解データ算出部とを有し、更に必要に応じてその他の部を有する。
【００１４】
＜ガウス・ザイデル法＞
図２に示すように、ガウス・ザイデル法は、連立方程式Ａｘ＝ｂについて、漸化式を反復計算することにより解ｘを求める方法である。漸化式では、更新済の解ベクトルｘの要素値を順次反映する。
ガウス・ザイデル法は、行列Ａの計算方法としては、行列Ａの右上部分「Ｕ」、行列Ａの左下部分「Ｌ」、行列Ａの対角要素部分「Ｄ」に分けて、それぞれを解ベクトルｘで掛け算を行う。
ガウス・ザイデル法は、行列Ａの「Ｕ」部分については、更新前のｘｉ（ｋ）を掛けるので、並列実行可能である。
ガウス・ザイデル法は、行列Ａの「Ｌ」部分については、更新後のｘｉ（ｋ＋１）を掛ける。即ち、ガウス・ザイデル法は、更新後のｘｉ（ｋ＋１）を参照し、行i−１の計算完了後に計算する。このように、ガウス・ザイデル法は、計算順序の依存性を有するため、並列実行が困難になる。
【００１５】
図３は、３×３行列のガウス・ザイデル法での計算の一例を示す図である。
ガウス・ザイデル法は、連立方程式の漸化式において、まず、ガウス・ザイデル法は、ｘ１（ｋ＋１）のステップを計算する。次に、ガウス・ザイデル法は、ｘ２（ｋ＋１）のステップを計算するときに、矢印部分では更新後のｘ１を参照して計算しなければならない。ガウス・ザイデル法は、次の、ｘ３（ｋ＋１）のステップを計算するときには、矢印部分で更新後のｘ１、ｘ２を参照して計算しなければならず、順番に処理しなければならない。このように、ガウス・ザイデル法は、連立方程式の左下三角形部分の非零要素の位置により、計算順序が決まってしまう。
【００１６】
＜逐次加速緩和法（ＳＯＲ法）＞
ＳＯＲ法は、ｎ元連立一次方程式Ａｘ＝ｂを反復法で解く手法の一つであり、ガウス・ザイデル法に加速パラメータωを導入してその修正量を拡大することで、更に高速化を図った手法である。
加速パラメータωは、下記数式で与えられる。
【数１】
2018206078.tif 000003
ＳＯＲ法は、収束するためには、０＜ω＜２であることが必要である。１＜ω＜２の場合は加速緩和法（ＳＯＲ法）、０＜ω＜１の場合は減速緩和法、ω＝１の場合はガウス・ザイデル法に対応する。
【００１７】
ここで、図４は、従来の並列処理装置全体の一例を示すブロック図である。この図４の並列処理装置２００は、右側並列演算結果データ算出部２０１、及び左側逐次演算結果データ算出部２０２を備える。
【００１８】
次に、図５は、従来の並列処理装置全体の処理の流れの一例を示すフローチャートである。以下、図４を参照して、従来の並列処理装置全体の処理の流れについて説明する。
【００１９】
ステップＳ１１では、右側並列演算結果データ算出部２０１は、所定行データの行列データの対角要素の右側を各行毎に並列演算して右側並列演算結果データを算出すると、処理をＳ１２に移行する。
【００２０】
ステップＳ１２では、左側逐次演算結果データ算出部２０２は、所定行データの行列データの対角要素の左側を各行毎に逐次演算して左側逐次演算結果データを算出すると、本処理を終了する。
【００２１】
図６は、従来のガウス・ザイデル法の一例を説明するための説明図である。
図６中Ｌ、Ｕは、行列サイズＮの行列ＡにおけるブロックサイズＮｂの部分行列である。
Ｄは、行列サイズＮの行列ＡにおけるブロックサイズＮｂの対角行列である。
ｘｎ、ｘｏは、解ベクトルｘの部分解ベクトルである。
ブロックサイズＮｂは、固定ではなく、行列Ａ中には零要素と、零以外の非零要素とがあり、計算するときは非零要素により計算するので、一律に決めることは困難であり、非零要素の数に応じて適宜設定される。
【００２２】
図７は、従来のガウス・ザイデル法による処理の流れの一例を示すフローチャートである。以下、図４及び図６を参照して、従来のガウス・ザイデル法による処理の流れについて説明する。
【００２３】
ステップＳ２１では、右側並列演算結果データ算出部２０１が、部分行列Ｕ、Ｌ、及び部分解ベクトルｘｎ、ｘｏを選択すると、処理をＳ２２に移行する。
【００２４】
ステップＳ２２では、右側並列演算結果データ算出部２０１は、部分行列Ｕと部分解ベクトルｘｏの行列ベクトル積を並列計算し、演算結果を一時ベクトルｔｍｐに格納すると、処理をＳ２３に移行する。
【００２５】
ステップＳ２３では、左側逐次演算結果データ算出部２０２は、部分行列Ｌと算出済み部分解ベクトルｘｎの行列ベクトル積を逐次計算し、解ベクトルｘｎを更新すると、処理をＳ２４に移行する。
【００２６】
ステップＳ２４では、左側逐次演算結果データ算出部２０２は、全ての解ベクトルｘを更新したか否かを判定する。左側逐次演算結果データ算出部２０２は、全ての解ベクトルｘを更新していないと判定すると、次の部分行列Ｕ、Ｌ、及び部分解ベクトルｘｎ、ｘｏを選択し、処理をＳ２２に移行する。一方、左側逐次演算結果データ算出部２０２は、全ての解ベクトルｘを更新したと判定すると、本処理を終了する。
【００２７】
図８は、従来の１６×１６行列のガウス・ザイデル法での計算の一例を示す図である。
並列実行部は、部分行列Ｕを並列計算する。
逐次実行部は、部分行列Ｌを逐次計算し、ｘ’９を計算した後、ｘ’１０を計算する時には、ｘ’９を使用して、ｘ’１０を算出する。ｘ’１１を計算する時には、ｘ’９、ｘ’１０を使用してｘ’１１を算出する。したがって、逐次実行部は、矢印の部分を順番に実施するという制約を有する。
【００２８】
まず、（１）の処理では、並列実行部は、部分行列Ｕの行列ベクトル積を並列計算し、その結果を一時ベクトルｔｍｐに格納する。
次に、（２）の処理では、逐次実行部は、部分行列Ｌの行列ベクトル積、及び一時ベクトルｔｍｐの結果から、解ベクトルｘの更新を逐次計算する。そのため、逐次実行部は、演算が進むにつれ、逐次計算される部分行列Ｌが大きくため、逐次計算の計算量が増加する。逐次実行部は、逐次計算の計算量が増えることで並列度が下がり、高速に処理を実行することが困難になる。
【実施例】
【００２９】
以下に、本発明の実施例について図面を用いて具体的に説明するが、本発明は、この実施例に何ら限定されるものではない。
【００３０】
（実施例１）
図９は、並列処理装置のハードウェア構成の一例を示す図である。この図９の並列処理装置１００の後述する主記憶装置及び記憶装置には、並列演算プログラムが記録されており、後述の演算装置が当該プログラムを読出して実行することにより、後述の右側演算結果データ算出部１１１、左側演算結果データ算出部１１２、及び対角要素解データ算出部１１３として動作する。
【００３１】
並列処理装置１００は、共有メモリ型計算機であり、相互接続ネットワーク１０４により互いに接続される演算装置１  １０１、演算装置２  １０２、・・・、演算装置ｎ  １０３、主記憶装置１０５、及び記憶装置１０６を備え、入力部１０７及び出力部１０８が接続されている。
【００３２】
演算装置１  １０１、演算装置２  １０２、・・・、演算装置ｎ  １０３は、主記憶装置１０５や記憶装置１０６に格納された右側演算結果データ算出部１１１、左側演算結果データ算出部１１２、及び対角要素解データ算出部１１３の各種プログラムを実行するユニットである。
【００３３】
相互接続ネットワーク１０４は、複数の機器がケーブルなどを介して接続され、これらの間でデータの中継及び転送を行い、機器間で通信ができるようにする装置である。
【００３４】
主記憶装置１０５は、例えば、ＲＯＭ（Ｒｅａｄ  Ｏｎｌｙ  Ｍｅｍｏｒｙ）、ＲＡＭ（Ｒａｎｄｏｍ  Ａｃｃｅｓｓ  Ｍｅｍｏｒｙ）等を含み、並列処理装置１００を構成する各処理のプログラムとデータが格納される。
【００３５】
記憶装置１０６としては、例えば、磁気ディスク装置、光ディスク装置、光磁気ディスク装置などが挙げられる。この記憶装置１０６に上述の各処理のプログラムとデータを保存しておき、必要に応じて、それらを主記憶装置１０５にロードして使用することもできる。
【００３６】
入力部１０７としては、例えば、キーボード、マウス、ポインティングデバイス、タッチパネルなどが挙げられ、利用者からの指示の入力に用いられる。また、可搬記録媒体を駆動してその記録内容の入力に用いられる。
【００３７】
出力部１０８としては、例えば、ディスプレイやプリンタ等であり、並列処理装置１００の利用者への処理結果等の表示に用いられる。
【００３８】
なお、図９には示していないが、演算処理の高速化のために、ＧＰＵ（Ｇｒａｐｈｉｃｓ  Ｐｒｏｃｅｓｓｉｎｇ  Ｕｎｉｔ）やＦＰＧＡ（Ｆｉｅｌｄ−Ｐｒｏｇｒａｍｍａｂｌｅ  Ｇａｔｅ  Ａｒｒａｙ）などのアクセラレータを利用できる構成としてもよい。
【００３９】
次に、図１０は、本発明の並列処理装置全体の一例を示すブロック図である。この図１０の並列処理装置１００は、右側演算結果データ算出部１１１、左側演算結果データ算出部１１２、及び対角要素解データ算出部１１３を備えている。
右側演算結果データ算出部１１１は、所定行データの行列データの対角要素の右側を各行毎に並列演算して右側演算結果データを算出する。
左側演算結果データ算出部１１２は、所定行データの行列データの対角要素の左側を算出済み解データを用い各行毎に並列演算して左側演算結果データを算出する。
対角要素解データ算出部１１３は、右側演算結果データと左側演算結果データとを用いて、所定行データの行列データの対角要素の解データを算出する。
【００４０】
ここで、右側演算結果データ算出部１１１、左側演算結果データ算出部１１２、及び対角要素解データ算出部１１３の構成は、本発明の並列処理装置に該当する。右側演算結果データ算出部１１１、左側演算結果データ算出部１１２、及び対角要素解データ算出部１１３を実施する処理は、本発明の並列演算方法に該当する。右側演算結果データ算出部１１１、左側演算結果データ算出部１１２、及び対角要素解データ算出部１１３の処理をコンピュータに実行させるプログラムは、本発明の並列演算プログラムに該当する。
【００４１】
ここで、図１１は、本発明の並列処理装置全体の処理の流れの一例を示すフローチャートである。以下、図１０を参照して、本発明の並列処理装置全体の処理の流れについて説明する。
【００４２】
ステップＳ３１では、右側演算結果データ算出部１１１は、所定行データの行列データの対角要素の右側を各行毎に並列演算して右側演算結果データを算出すると、処理をＳ３２に移行する。
【００４３】
ステップＳ３２では、左側演算結果データ算出部１１２は、所定行データの行列データの対角要素の左側と算出済み解データを用い各行毎に並列演算して左側演算結果データを算出すると、処理をＳ３３に移行する。
【００４４】
ステップＳ３３では、対角要素解データ算出部１１３は、右側演算結果データと左側演算結果データとを用いて、所定行データの行列データの対角要素の解データを算出すると、本処理を終了する。
【００４５】
図１２は、本発明のガウス・ザイデル法の一例を説明するための説明図である。
まず、ガウス・ザイデル法においては、並列処理装置は、漸化式のＬ部分をＬｐ部分とＬｓ部分とに分割する。
次に、並列処理装置は、Ｕ部分、及びＬｐ部分を並列計算する。
次に、並列処理装置は、Ｌｓ部分、Ｄ部分、及び解ベクトルｘの更新を逐次計算する。
【００４６】
図１２中Ｌｐ、Ｌｓ、及びＵは、行列サイズＮの行列ＡにおけるブロックサイズＮｂの部分行列である。
Ｄは、行列サイズＮの行列ＡにおけるブロックサイズＮｂの対角要素である。
ブロックサイズＮｂは、固定ではなく、行列Ａ中には零要素と、零以外の非零要素とがあり、計算するときは非零要素により計算するので、一律に決めることは困難であり、非零要素の数に応じて適宜設定される。
部分行列Ｕは、所定行データの行列データの対角要素Ｄの右側である。
部分行列Ｌｐは、所定行データの行列データの対角要素Ｄの左側である。
解ベクトルｘは、部分解ベクトルｘｎ１、ｘｏ１、ｘｏ２からなる。部分解ベクトルｘｎ１は、算出済み解データを表す。部分解ベクトルｘｏ１は、所定行データの対角要素Ｄの解データである。
【００４７】
図１３は、本発明のガウス・ザイデル法による処理の流れの一例を示すフローチャートである。以下、図１０及び図１２を参照して、本発明のガウス・ザイデル法による処理の流れについて説明する。
【００４８】
ステップＳ４１では、右側演算結果データ算出部１１１が、部分行列Ｕ、Ｌｐ、Ｌｓ、及び部分解ベクトルｘｎ１、ｘｏ１、ｘｏ２を選択すると、処理をＳ４２に移行する。
【００４９】
ステップＳ４２では、右側演算結果データ算出部１１１は、部分行列Ｕと部分解ベクトルｘｏ２の行列ベクトル積を並列計算し、演算結果を一時ベクトルｔｍｐに格納すると、処理をＳ４３に移行する。
【００５０】
ステップＳ４３では、左側演算結果データ算出部１１２は、部分行列Ｌｐと部分解ベクトルｘｎ１の行列ベクトル積を並列計算し、演算結果を一時ベクトルｔｍｐに加算すると、処理をＳ４４に移行する。
【００５１】
ステップＳ４４では、対角要素解データ算出部１１３は、部分行列Ｌｓと部分解ベクトルｘｏ１の行列ベクトル積を逐次計算し、演算結果を一時ベクトルｔｍｐに加算して、部分解ベクトルｘｏ１を更新すると、処理をＳ４５に移行する。
【００５２】
ステップＳ４５では、対角要素解データ算出部１１３は、全ての解ベクトルｘを更新したか否かを判定する。対角要素解データ算出部１１３は、全ての解ベクトルｘを更新していないと判定すると、次の部分行列Ｕ、Ｌｐ、Ｌｓ、及び部分解ベクトルｘｎ１、ｘｏ１、ｘｏ２を選択し、処理をＳ４２に移行する。一方、対角要素解データ算出部１１３は、全ての解ベクトルｘを更新したと判定すると、本処理を終了する。
【００５３】
図１４は、本発明の１６×１６行列のガウス・ザイデル法での計算の一例を示す図である。
並列実行部は、部分行列Ｕを並列計算する。部分行列Ｌｐを並列計算する。すると、ｘ’１〜ｘ’８は、ｘ’９を計算する前に計算されている。
逐次実行部は、部分行列Ｌｓを逐次計算する。
【００５４】
まず、（１）の処理では、ガウス・ザイデル法での並列処理装置は、部分行列Ｕの行列ベクトル積を並列計算し、その結果を一時ベクトルｔｍｐに格納する。
次に、（２）の処理では、並列処理装置は、部分行列Ｌｐの行列ベクトル積を並列計算し、その結果を一時ベクトルｔｍｐに加算する。
次に、（３）の処理では、並列処理装置は、部分行列Ｌｓ、行列ベクトル積、及び一時ベクトルｔｍｐの結果から、解ベクトルｘの更新を逐次計算する。
図１４は、図８と比べて、ＬがＬｐとＬｓに分けられ、Ｌｐが並列計算でき、Ｌｓのみが逐次計算となるで、並列化により高速に実行することができる。
【００５５】
ここで、再び、図１に示すガウス・ザイデル法による処理について説明する。
まず、ガウス・ザイデル法は、ｉ＝ｉ１行目の実行を開始する時点を考えると、ｉ＝ｉ１−１行目までの計算により、ベクトルｘの要素ｘ［０］〜ｘ［ｉ１−１］の値は更新済みである。以下、既に値が更新されていることを、値は「ＮＥＷ」であると表現する。
一方、要素ｘ［ｉ１］〜ｘ［ｎ−１］はまだ値が更新されていない。以下、まだ値が更新されていないことを、値は「ＯＬＤ」であると表現する。
【００５６】
次に、ガウス・ザイデル法は、ｉ＝ｉ１行目を実行する場合、図１中の（ｄ）のｘ［Ａｃｏｌ［ｉ］［ｊ］］は、０≦Ａｃｏｌ［ｉ１］［ｊ］≦ｉ１−１ならば、ｘ［Ａｃｏｌ［ｉ］［ｊ］］の値は「ＮＥＷ」である。
一方、ガウス・ザイデル法は、ｉ１＋１≦Ａｃｏｌ［ｉ１］［ｊ］≦ｎ−１ならば、図１中の（ｄ）のｘ［Ａｃｏｌ［ｉ］［ｊ］］の値は「ＯＬＤ」である。即ち、ガウス・ザイデル法は、ｉ＝ｉ１行目の図１中の（ｄ）の積差計算をｘ［Ａｃｏｌ［ｉ］［ｊ］］が「ＮＥＷ」だけである部分、ｘ［Ａｃｏｌ［ｉ］［ｊ］］が「ＯＬＤ」だけである部分の２つに分けることができる。
「ＯＬＤ」だけである部分の計算は、要素ｘ［ｉ１＋１］〜ｘ［ｎ−１］が「ＯＬＤ」である間、即ち、要素ｘ［ｉ＋１］が「ＮＥＷ」となる前であれば、いつでも（例えば、要素ｘ［０］〜ｘ［ｎ−１］がまだ「ＯＬＤ」である時から）計算可能である。
【００５７】
以上を踏まえて、疎行列Ａを図１５に示すように分割する。この図１５は、疎行列Ａを行方向にブロック０〜３の４つのブロックに分割する。更に、各ブロックを、対角線Ｄに沿った下三角形領域Ｔ０、Ｔ１、Ｔ２、Ｔ３と、それ以外の右側領域Ｒ０、Ｒ１、Ｒ２、Ｒ３、及び左側領域Ｌ１、Ｌ２、Ｌ３に分割する。
【００５８】
図１５の疎行列Ａにおいて、ブロック０、ブロック１、ブロック２、及びブロック３の順に計算を実行する。
まず、ブロック０において、疎行列Ａの下三角形領域Ｔ０の非零要素との計算に使用される要素ｘ［ｉ］は「ＮＥＷ」である。
疎行列Ａの下三角形領域Ｔ０の右側領域Ｒ０の非零要素との計算に使用される要素ｘ［ｉ］は「ＯＬＤ」であるため、各行の計算を「ＮＥＷ」と「ＯＬＤ」に分ければ、ブロック０内の実行順序を下三角形領域の右側領域Ｒ０、下三角形領域Ｔ０の順とすることができ、かつ、下三角形領域の右側領域Ｒ０は全行を同時に並列実行可能である。一方、下三角形領域Ｔ０は逐次実行可能である。
【００５９】
次に、ブロック１について、下三角形領域Ｔ１の右側領域Ｒ１は、「ＯＬＤ」であるため、全行を同時に並列実行可能である。
下三角形領域の左側領域Ｌ１は、既にブロック０で「ＮＥＷ」となっているベクトルｘの要素を参照できるので、ブロック１内では、全行を同時に並列実行可能である。即ち、下三角形領域の右側領域Ｒ１及び左側領域Ｌ１は、ブロック１内では全行同時に並列実行可能であり、下三角形領域Ｔ１は逐次実行可能である。
その結果、ブロック１では、ブロック０と同様に、下三角形領域の右側領域Ｒ１及び左側領域Ｌ１、下三角形領域Ｔ１の順で実行し、かつ下三角形領域の右側領域Ｒ１及び左側領域Ｌ１は全行を並列実行可能である。
【００６０】
次に、ブロック２について、下三角形領域Ｔ２の右側領域Ｒ２は、「ＯＬＤ」であるため、全行を同時に並列実行可能である。
下三角形領域Ｔ２の左側領域Ｌ２は、既にブロック１で「ＮＥＷ」となっているベクトルｘの要素を参照できるので、ブロック２内では、全行を同時に並列実行可能である。即ち、下三角形領域の右側領域Ｒ２及び左側領域Ｌ２は、ブロック２内では全行同時に並列実行可能であり、下三角形領域Ｔ２は逐次実行可能である。
その結果、ブロック２では、ブロック０、１と同様に、下三角形領域の右側領域Ｒ２及び左側領域Ｌ２、下三角形領域Ｔ２の順で実行し、かつ下三角形領域の右側領域Ｒ２及び左側領域Ｌ２は全行を並列実行可能である。
【００６１】
次に、ブロック３について、下三角形領域Ｔ３の右側領域Ｒ３は、「ＯＬＤ」であるため、全行を同時に並列実行可能である。
下三角形領域Ｔ３の左側領域Ｌ３は、既にブロック２で「ＮＥＷ」となっているベクトルｘの要素を参照できるので、ブロック３内では、全行を同時に並列実行可能である。即ち、下三角形領域の右側領域Ｒ３及び左側領域Ｌ３は、ブロック３内では全行同時に並列実行可能であり、下三角形領域Ｔ３は逐次実行可能である。
その結果、ブロック３では、ブロック０、１、２と同様に、下三角形領域の右側領域Ｒ３及び左側領域Ｌ３、下三角形領域Ｔ３の順で実行し、かつ下三角形領域の右側領域Ｒ３及び左側領域Ｌ３は全行を並列実行可能である。
【００６２】
したがって、ベクトルｘの更新順序（要素ｘ［０］からｘ［ｎ−１］まで順に更新していく）を変更せずに、図１６に示す（１）〜（８）の順に、下三角形領域Ｔ０、Ｔ１、Ｔ２、Ｔ３を除いた右側領域Ｒ０、Ｒ１、Ｒ２、Ｒ３及び左側領域Ｌ１、Ｌ２、Ｌ３は、並列実行可能である。
【００６３】
具体的に、プロセッサ０〜Ｎ−１のＮ個のプロセッサは、図１６に示すような疎行列Ａを並列実行する。なお、図１６は、Ｎ＝４の場合である。
まず、（１）プロセッサ０〜Ｎ−１は、下三角形領域Ｔ０の右側領域Ｒ０を並列実行する。
次に、（２）プロセッサ０は、下三角形領域Ｔ０を逐次実行する。
次に、（３）プロセッサ０〜Ｎ−１は、下三角形領域Ｔ１の右側領域Ｒ１及び左側領域Ｌ１を並列実行する。
次に、（４）プロセッサ０は、下三角形領域Ｔ１を逐次実行する。
次に、（５）プロセッサ０〜Ｎ−１は、下三角形領域Ｔ２の右側領域Ｒ２及び左側領域Ｌ２を並列実行する。
次に、（６）プロセッサ０は、下三角形領域Ｔ２を逐次実行する。
次に、（７）プロセッサ０〜Ｎ−１は、下三角形領域Ｔ３の右側領域Ｒ３及び左側領域Ｌ３を並列実行する。
最後に、（８）プロセッサ０は、下三角形領域Ｔ３を逐次実行する。
以上説明したように、実施例１によれば、疎行列に対するガウス・ザイデル法を並列化して高速に実行することができる。なお、本実施例では、プロセッサの並列処理について説明したが、ノードの並列処理、プロセッサ内のコアの並列処理についても適用できる。
【００６４】
（実施例２）
実施例１は、図１６に示すように、下三角形領域Ｔ０、Ｔ１、Ｔ２、及びＴ３をプロセッサ０のみで実行し、残りのプロセッサ１〜Ｎ−１はアイドル状態となっている。そこで、実施例２では、プロセッサ０〜Ｎ−１のＮ個のプロセッサは、図１７に示すように、疎行列Ａを並列実行する。なお、図１７は、Ｎ＝４の場合である。
まず、（１）プロセッサ０〜Ｎ−１は、下三角形領域Ｔ０の右側領域Ｒ０を並列実行する。
次に、（２）プロセッサ０は、下三角形領域Ｔ０を逐次実行すると同時に、プロセッサ１〜Ｎ−１は、下三角形領域Ｔ１の右側領域Ｒ１を並列実行する。
次に、（３）プロセッサ０〜Ｎ−１は、下三角形領域Ｔ１の左側領域Ｌ１を並列実行する。
次に、（４）プロセッサ０は、下三角形領域Ｔ１を逐次実行すると同時に、プロセッサ１〜Ｎ−１は、下三角形領域Ｔ２の右側領域Ｒ２を並列実行する。
次に、（５）プロセッサ０〜Ｎ−１は、下三角形領域Ｔ２の左側領域Ｌ２を並列実行する。
次に、（６）プロセッサ０は、下三角形領域Ｔ２を逐次実行すると同時に、プロセッサ１〜Ｎ−１は、下三角形領域Ｔ３の右側領域Ｒ３を並列実行する。
次に、（７）プロセッサ０〜Ｎ−１は、下三角形領域Ｔ３の左側領域Ｌ３を並列実行する。
最後に、（８）プロセッサ０は、下三角形領域Ｔ３を逐次実行する。
このようにすることにより、プロセッサ０のみで実行するのは下三角形領域Ｔ３のみとなる。
実施例２によれば、アイドル状態のプロセッサをより少なくでき、疎行列に対するガウス・ザイデル法を並列化して高速に実行することができる。
【００６５】
（実施例３）
実施例１において、ガウス・ザイデル法の代わりにＳＯＲ法を用いた以外は、実施例１と同様にして、プロセッサ０〜Ｎ−１のＮ個のプロセッサは、図１６に示すように、疎行列Ａを並列実行する。なお、図１６は、Ｎ＝４の場合である。
ＳＯＲ法は、実施例１のガウス・ザイデル法に加速パラメータω（１＜ω＜２）を導入した点が異なるだけである。
実施例３によれば、アイドル状態のプロセッサを少なくでき、疎行列に対するＳＯＲ法を並列化してより高速に実行することができる。
【００６６】
（実施例４）
実施例２において、ガウス・ザイデル法の代わりにＳＯＲ法を用いた以外は、実施例２と同様にして、図１７に示すように、疎行列Ａをプロセッサ０〜Ｎ−１のＮ個のプロセッサで並列計算を実行した。なお、図１７は、Ｎ＝４の場合である。
ＳＯＲ法は、実施例１のガウス・ザイデル法に加速パラメータω（１＜ω＜２）を導入した点が異なるだけである。
実施例４によれば、疎行列に対するＳＯＲ法を並列化してより高速に実行することができる。
【００６７】
以上の実施形態に関し、更に以下の付記を開示する。
（付記１）
ガウス・ザイデル法及び逐次加速緩和法（ＳＯＲ法）の少なくともいずれかで行列データを用いて行列連立一次方程式を解く際に、前記行列データの所定行データずつ演算して各行毎に解データを出力する並列処理装置において、
前記所定行データの行列データの対角要素の右側を各行毎に並列演算して右側演算結果データを算出する右側演算結果データ算出部と、
前記所定行データの行列データの対角要素の左側と算出済み解データを用い各行毎に並列演算して左側演算結果データを算出する左側演算結果データ算出部と、
前記右側演算結果データと前記左側演算結果データとを用いて、前記所定行データの行列データの対角要素の解データを算出する対角要素解データ算出部と、
を有することを特徴とする並列処理装置。
（付記２）
前記対角要素の解データを算出する際に、次の所定行データの行列データの対角要素の右側を各行毎に並列演算して右側演算結果データを算出する付記１に記載の並列処理装置。
（付記３）
ガウス・ザイデル法及び逐次加速緩和法（ＳＯＲ法）の少なくともいずれかで行列データを用いて行列連立一次方程式を解く際に、前記行列データの所定行データずつ演算して各行毎に解データを出力する並列演算方法において、
前記所定行データの行列データの対角要素の右側を各行毎に並列演算して右側演算結果データを算出し、
前記所定行データの行列データの対角要素の左側と算出済み解データを用い各行毎に並列演算して左側演算結果データを算出し、
前記右側演算結果データと前記左側演算結果データとを用いて、前記所定行データの行列データの対角要素の解データを算出することを特徴とする並列演算方法。
（付記４）
前記対角要素の解データを算出する際に、次の所定行データの行列データの対角要素の右側を各行毎に並列演算して右側演算結果データを算出する付記３に記載の並列演算方法。
（付記５）
ガウス・ザイデル法及び逐次加速緩和法（ＳＯＲ法）の少なくともいずれかで行列データを用いて行列連立一次方程式を解く際に、前記行列データの所定行データずつ演算して各行毎に解データを出力する並列演算プログラムにおいて、
前記所定行データの行列データの対角要素の右側を各行毎に並列演算して右側演算結果データを算出し、
前記所定行データの行列データの対角要素の左側と算出済み解データを用い各行毎に並列演算して左側演算結果データを算出し、
前記右側演算結果データと前記左側演算結果データとを用いて、前記所定行データの行列データの対角要素の解データを算出する
処理をコンピュータに実行させることを特徴とする並列演算プログラム。
（付記６）
前記対角要素の解データを算出する際に、次の所定行データの行列データの対角要素の右側を各行毎に並列演算して右側演算結果データを算出する付記５に記載の並列演算プログラム。
【符号の説明】
【００６８】
１００      並列処理装置
１０１      演算装置１
１０２      演算装置２
１０３      演算装置ｎ
１０４      相互接続ネットワーク
１０５      主記憶装置
１０６      記憶装置
１０７      入力部
１０８      出力部
１１１      右側演算結果データ算出部
１１２      左側演算結果データ算出部
１１３      対角要素解データ算出部
【図１】
2018206078.tif 000004
【図２】
2018206078.tif 000005
【図３】
2018206078.tif 000006
【図４】
2018206078.tif 000007
【図５】
2018206078.tif 000008
【図６】
2018206078.tif 000009
【図７】
2018206078.tif 000010
【図８】
2018206078.tif 000011
【図９】
2018206078.tif 000012
【図１０】
2018206078.tif 000013
【図１１】
2018206078.tif 000014
【図１２】
2018206078.tif 000015
【図１３】
2018206078.tif 000016
【図１４】
2018206078.tif 000017
【図１５】
2018206078.tif 000018
【図１６】
2018206078.tif 000019
【図１７】
2018206078.tif 000020