JP2006085619A - 帯係数行列を持つ連立１次方程式の解法プログラム - Google Patents

Abstract

【課題】 圧縮モードで格納された帯行列を用いて、並列処理によってＬＵ分解を効率化し、連立１次方程式の計算を高速にする。
【解決手段】 解法プログラムは、複数の列から構成される列ブロックを並列的にＬＵ分解して作業領域に格納する手順１、１の結果に対して各列の左側での行の入替をキャンセルして圧縮モードの配列にコピーバックする手順２、１の結果に対応して帯行列の更新によって壊れる可能性のある部分を退避する手順３、１の結果を用いて帯行列を並列的に更新する手順４、４の結果上に３で退避した部分を戻す手順５を計算機に実行させる。
【選択図】図１

Description

本発明は係数行列として疎行列、すなわち０でない要素が少ない行列の代表例としての帯行列を持つ連立１次方程式の解法に係り、さらに詳しくはそのような連立１次方程式を共用メモリ型スカラ並列計算機によって解くための連立１次方程式の解法プログラムに関する。

連立１次方程式を計算機によって解く場合には、連立１次方程式を行列によって表示し、行列のＬＵ分解などの演算を行って解を求めやすい形式に変形し、方程式の解を求めるガウスの消去法を基本とする方法が用いられる。

すなわち、連立１次方程式は係数を表わす行列と変数を表わす列ベクトルの積が定数列ベクトルと等しくなるという形式に記述することができる。ここでＬＵ分解を行って、係数を表わす行列を上三角行列と下三角行列とに分解し、前進代入（前進消去）と後退代入という方法を用いることによって、連立１次方程式の解を求めることができる。したがって連立１次方程式を解くためには係数行列をＬＵ分解することが重要な処理となる。このＬＵ分解を共用メモリ型スカラ並列計算機を用いて効率的に並列処理する従来技術として、出願人の次の特許文献がある。
特開２００２−１６３２４６号公報 「共有メモリ型スカラ並列計算機における並列行列処理方法、及び記録媒体」

この文献には、ＬＵ分解すべき行列のうちで左側の複数の列に対応する左上部の対角部分のブロックＤと、Ｄの下側にある列方向のブロックとを取り、下側にある列方向のブロックを、例えばＬ１からＬ３の３つに分け、３つのプロセッサのそれぞれにＤ＋Ｌ１、Ｄ＋Ｌ２、Ｄ＋Ｌ３を割り振り、並列にＬＵ分解演算を行い、その後対角部分のブロックＤの右側の複数の行によって構成されるブロックＵを更新し、さらにＬ１からＬ３とＵとを用いて行列の残りの部分を更新する処理を繰り返すことによって、処理の効率化を実現できる演算方法が開示されている。

係数行列の中で０でない要素の数が少ない疎行列のうちで、特に主対角線の周囲にしか０でない要素が存在しない帯行列（バンド行列）を係数行列として持つ連立１次方程式の解法において、従来ではガウスの消去法をベースとしたＬＵ分解を行う方法が用いられている。このような従来の解法では行列の要素をメモリに格納する場合に、０でない要素が存在する帯の部分だけを効率的に格納するために、帯の部分以外の０の要素の格納を省略する圧縮格納モードが用いられる。またＬＵ分解の解の安定性を高めるために部分ピポットが採用されるが、圧縮格納モードでの格納領域を小さくするために、各列でピポットを用いて行の入れ替えを行う場合に列の右側のみで行の入れ替えを行う方法が用いられていた。さらにＬＵ分解における更新の処理では、ベクトルの外積の形式の演算が用いられていた。

このような帯係数行列を持つ連立１次方程式の解法の従来技術は次の文献に記載されている。
Ｇ．Ｈ．Ｇｏｒｕｂ，Ｃ．Ｆ．Ｖａｎ Ｌｏａｎ：Ｍａｔｒｉｘ Ｃｏｍｐｕｔａｔｉｏｎｓ，３ｒｄ Ｅｄ．Ｔｈｅ Ｊｏｈｎｓ Ｈｏｐｋｉｎｓ Ｕｎｉｖｅｒｓｉｔｙ Ｐｒｅｓｓ，Ｂａｌｔｉｍｏｒｅ ａｎｄ Ｌｏｎｄｏｎ （１９９６）

一般にスカラ計算機は、ＣＰＵの演算性能は高いが、メモリへのアクセス性能が低い。このため、メモリアクセスの性能に依存するベクトルの外積をベースとする演算は性能が低く、ベクトル計算機による処理に比べて効率が悪いという第１の問題点があった。

第２に帯行列が０でない要素のみを格納する圧縮モードでメモリに格納されているため、例えば各列が行方向に１要素ずつずれて格納され、そのままの形式では行列積を用いた更新処理を行うことができず、また行列積を用いて更新しようとする場合にも入れ替えの必要な行における要素の数が圧縮モードの格納領域を越えた時に、その外側の要素の値が壊れてしまう可能性があるという問題点があった。

第３に前述のように部分ピポットの形式で列の右側のみで行の入れ替えを行うために、ブロック化した行列演算の形を利用する更新を行うことができないという問題点があった。

第４にＬＵ分解が終了して前進消去の処理を行う場合に、部分ピポットの形式で列の右側部分のみ行の入れ替えを行っているために、解ベクトルの入れ替えも行う必要があり、ベクトルとスカラの積の演算によって更新処理を行うことが必要となり、その演算を並列に処理する場合には並列処理そのもののオーバヘッドが大きくなり、並列処理の効率が悪化するという問題点があった。

本発明の課題は上述の問題点に鑑み、帯係数行列を持つ連立１次方程式の解法において、帯行列の圧縮格納モードによって生ずる問題点を解決して、ブロック化した行列演算の形式を利用することによってＬＵ分解の高速化を行うことと、ＬＵ分解結果に対する前進消去処理における並列化の効率を向上させることである。

図１は本発明の連立１次方程式解法プログラムの原理的な機能ブロック図である。同図は、帯係数行列を持つ連立１次方程式の解法プログラムであり、例えば共用メモリ型スカラ並列計算機によって使用されるものである。

図１において、プログラムでは１で帯行列内部の複数の列から構成される列ブロックが並列的にＬＵ分解されてその結果が作業領域に格納され、２でそのＬＵ分解結果に対して各列のピボット選択の結果として左側で行われていた行の入れ替えがキャンセルされ、そのキャンセル結果のデータが圧縮形式の帯行列にコピーバックされ、３でこの列ブロックのＬＵ分解結果に対応して入れ替えの必要な行の長さの最大値と圧縮形式の帯行列の配列領域の大きさとの関係に基づいて、行列の更新によって壊れる可能性のある部分が退避される。

その後帯行列の残りの部分の更新が行われる。すなわち前述の１で作業領域に格納された列ブロックのＬＵ分解結果を用いて、４で並列的に帯行列の更新が行われ、５ですでに３で退避された部分が帯行列の更新結果上に戻される。すなわち本発明のプログラムは１〜５の手順から成る行列演算を計算機に実行させることになる。

発明の実施の形態においては、この行列演算が帯行列のＬＵ分解のための演算であり、この演算では帯行列の最も左上部の前述の列ブロックに対応する行列演算の終了後に、その列ブロックの最上部の対角部分に含まれる行と列とを帯行列から除外し、除外後の行列の左上部の列ブロックに対応する行列演算を繰返し、最後に残った部分のＬＵ分解を行うこともできる。

また実施の形態においては、この残った部分のＬＵ分解の演算の終了後に、前述の帯行列の並列的更新において行われた行入れ替えの情報を利用して、複数の列から構成される列ブロックの各列の左側でピボット選択に対応する行の入れ替えを行うとともに、その左側での行の入れ替えに対応して連立方程式における定数ベクトルの要素の中でその列に対応する要素の入れ替えを行う手順と、その列ブロックの対角部分の下三角行列に対応する連立方程式を解く手順と、列ブロックの対角部分の下の行列を用いて行列ベクトル積の計算によってベクトルの更新を行う手順とを含む並列的行列演算をさらに計算機に実行させることもできる。

この場合、この並列的行列演算が帯行列のＬＵ分解結果を用いた前進代入の処理のための演算であり、帯行列の最も左上部の列ブロックに対応する並列的行列演算の終了後に、その列ブロックの最上部の対角部分に含まれる行と列とを帯行列から除外して、除外後の行列の左上部の列ブロックに対応する並列的行列演算を繰返し、最後に残った対角ブロックの下三角行列に対する連立方程式の解を求めることもできる。

さらに実施の形態においては、図１の４における帯行列の並列的更新手順において、キャンセル結果のデータがコピーバックされた圧縮形式の帯行列に対して、各列の要素の行の位置を補正して補正後の行列の更新を行うこともでき、またこの補正後の帯行列と列ブロックのＬＵ分解結果を用いて、列ブロック内の対角行列に対応する行ブロックを並列的に更新し、その行ブロックの更新結果を用いて列ブロックと行ブロックとに対応する行列を並列的に更新することもできる。

次に本発明の帯係数行列を持つ連立方程式の解法プログラムは、帯係数行列のＬＵ分解の演算の終了後に、帯係数行列の並列的更新において行われた行入れ替えの情報を利用して、複数の列から構成される列ブロックの各列の左側で、ピボット選択に対応する行の入れ替えを行うとともに、その行の入れ替えに対応して連立方程式における定数ベクトルの要素内でその列ブロックに対応する要素の入れ替えを行う手順と、列ブロックの対角部分の下三角行列に対応する連立方程式を解く手順と、列ブロックの対角部分の下の行列を用いて行列ベクトル積の計算によってベクトルの更新を行う手順とを含む行列演算を計算機に実行させるものである。

本発明によれば、帯行列を係数行列として持つ連立１次方程式の解法において、帯行列が圧縮格納モードでメモリに格納されている場合にも、係数行列のＬＵ分解を高速化することができ、また前進代入の処理においても並列処理の効率を向上させることができ、帯係数行列を持つ連立１次方程式の解を求める処理の高速化に寄与するところが大きい。

図２は、本発明のプログラムを用いて連立１次方程式の解を求めるための共用メモリ型スカラ並列計算機のハードウェア構成例を示すブロック図である。同図において共用メモリ型スカラ並列計算機を構成する複数のプロセッサ１０−１、１０−２、．．１０−ｎは、２次キャッシュメモリ１３−１、１３−２、．．１３−ｎを介して相互結合網１２に接続される。

各プロセッサ１０−１、１０−２、．．１０−ｎは、その内部、あるいは各２次キャッシュメモリと各プロセッサとの間に１次キャッシュメモリを備えている。また各プロセッサ１０−１、１０−２、．．１０−ｎによって共有されるメモリモジュール１１−１、１１−２、．．１１−ｎに対しては、相互結合網１２を介してプロセッサ１０−１、１０−２、．．１０−ｎがアクセスできるように接続されている。

プロセッサ１０−１、１０−２、．．１０−ｎがデータ処理を行う場合には、まずメモリモジュール１１−１、１１−２、．．１１−ｎから自プロセッサが担当するデータを２次キャッシュメモリ１３−１、１３−２、．．１３−ｎに格納し、さらにその中から処理単位となるデータを１次キャッシュメモリにコピーして処理を行うことになる。

処理が終ると、１次キャッシュメモリから２次キャッシュメモリに処理データが格納され、２次キャッシュメモリ内のデータに対する処理が全て終了すると、メモリモジュール１１−１、１１−２、．．１１−ｎの中で最初にデータを持ってきたメモリモジュールに対するデータの更新が行われる。このような処理を繰返すことによって複数のプロセッサによって並列処理が行われる。

この場合、各プロセッサが処理後のデータをメモリモジュールに書込み、次の処理のために再びメモリモジュールからデータを読み込む時にはプロセッサの間の同期をとる必要がある。すなわち、すべてのプロセッサがメモリモジュールのデータを更新し終わるまで他のプロセッサがメモリモジュールからデータを読み込まないようにする必要がある。このようなプロセッサ間の処理の同期をバリア同期という。

図３は、本発明で対象とする連立１次方程式の係数行列としての帯行列の説明図である。同図において帯行列はＮ行、Ｎ列の正方行列であり、その主対角線の周囲にのみ０でない要素が存在し、主対角線を中心とするある帯の幅の外側の要素はすべて０となっている疎行列である。例えば、中央付近の１つの列を取り、対角要素を除いた下側のその列内で０でない可能性のある要素の数を下バンド幅ｎｈ１、対角要素を除いた上側の列内の０でない要素の数を上バンド幅ｎｈ２と呼ぶことにする。

図４は、図３の帯行列の圧縮モードにおける格納形式例の説明図である。同図においておいて１番上側の下バンド幅ｎｈ１を行数として持つ領域は、後述する行の入れ替えによって０でない値が入る可能性のある領域であり、最初に帯行列を格納した段階では全ての要素が０となっている領域である。その下の行数ｎｈ２の領域は上バンド幅の要素を格納する領域であり、その下に対角要素を格納する１行の領域があり、さらにその下に下バンド幅ｎｈ１の要素を格納する領域が存在する。

図３の帯行列を図４のように圧縮モードで格納する場合には、例えば図３で最も左側の列から１列ずつ図４の一番上の行数ｎｈ１の領域を除いて格納していくことにする。すなわち最も左側の列は図３において対角要素から始まる列であり、図４の対角要素を格納する領域、すなわち行数１の領域から始めて下方向にｎｈ１の行数の領域に０でない（可能性のある）要素が格納される。２列目では、対角要素の上に１つの要素が存在するためにその要素が行数ｎｈ２の領域に１つだけ一番下に格納され、その後対角要素、およびその下の要素が次々と格納される。

格納される列が上のバンド（帯）境界に達した以後は、ｎｈ２の行数の領域、およびｎｈ１の行数の領域には列ごとに全て０でない（可能性のある）要素が格納され、格納される要素数はｎｈ１＋ｎｈ２＋１で一定となる。

さらに列の下端が図３の行列の最下行に達した以後は、図４の下側の行数ｎｈ１の領域に格納される要素は１つずつ減り、最も右側の列では格納される要素は上バンド側の要素と対角要素のみとなる。したがって図４において、行数ｎｈ２の領域の左上側、および下の行数ｎｈ１の領域の右下側の三角の領域には、不定の値、例えば全て０が格納されることになる。なお圧縮モードでの帯行列の格納形式としては、このような列単位での格納でなく、行単位での格納など、各種の格納形式を用いることができる。

図５は、ＬＵ分解を行うための列ブロックの作業領域への格納方法の説明図である。基本的にはこの列ブロックの格納では、最初に図３において最も左上側から対角要素を含む列ブロックの格納が行われるが、ここでは一般的に任意の位置の対角要素から始まる列ブロックの作業領域への格納について説明する。

本実施形態における列ブロックは、対角要素を最も左上の要素とする複数列の要素によって構成されるブロックであり、一般的には上側の対角部分と下側の対角要素を含まない部分とによって構成される。

この列ブロックを構成する複数の列内の要素は、図４、および図５に説明した圧縮モードでは、図３に示すような正しい列と行の位置関係ではなく、一列ごとに行が１つずつずれた形式で格納されている。すなわち図５において、作業領域Ｗにコピーすべき列ブロックにおける最初の列の要素は対角要素から始まり、その次の列の要素は対角要素の１つ上の行の要素から始まることになる。このように順次１つずつ行のずれた要素を、このずれを補正した形式で、作業領域Ｗとしての長方形の配列に行と列が正常な位置関係となるような形式で格納する。このように行と列の関係がずれのない正常な位置関係に補正された後に列ブロックのＬＵ分解が行われるために、そのＬＵ分解の結果を用いた行列の更新を行列積の計算によって行うことが可能となり、これによって並列処理による処理効率が向上する。

図６は、この列ブロックの帯係数行列上での位置の説明図である。同図において左上の角が主対角線上にある要素であり、複数の列で構成される長方形のブロックが、作業領域Ｗに格納される列ブロックを示す。この列ブロックを含む横長の実線の長方形の領域は、図５では下側の図で、左上の縦長の長方形とその右側の横長の長方形とに相当することになる。なお図５において斜めの実線、および点線は全て平行であり、その傾きは図４で説明した斜め方向の実線の傾きと一致する。また図６内の退避必要領域については後述する。

図７は、作業領域Ｗにコピーされた列ブロックの左下三角部分への０の設定の説明図である。この下三角領域は、図６で示すように下側のバンドの境界を越えた領域に相当し、この領域の要素は全て０であるために図７に示すように作業領域Ｗにコピーされた列ブロックの左下三角部分には０が設定される。

以上のようにして作業領域Ｗにコピーされた列ブロックに対するＬＵ分解が、前述の特許文献１において開示された方法を用いて並列的に実行される。この並列処理の詳細については特許文献１に記載されているが、その概要を図８を用いて説明する。

図８においては、例えば３つのプロセッサによる並列処理を行うものとし、列ブロックの対角部分をＤとし、その下の部分を１次元目、すなわち行数で均等に分割したＬ１、Ｌ２、およびＬ３を３つのスレッド（プロセッサ）Ｔ１からＴ３に割当て、各スレッドがＤ＋Ｌ１、Ｄ＋Ｌ２、Ｄ＋Ｌ３の演算を行う、すなわち対角部分を冗長に演算することによって列ブロックのＬＵ分解が行われる。このＬＵ分解において部分ピボットを用いた行の入れ替えでは列の右側だけでなく、ブロック幅全体に対して行の入れ替えが行われる。

ＬＵ分解された結果は圧縮モードの帯行列格納領域、すなわち図４で説明した領域内の対応する列にコピーバックされる。このコピーバックに先立って、ＬＵ分解処理の中で行われた行の入れ替えのうちで各列の左側における行の入れ替えがキャンセルされてから、コピーバックが行われる。

図９は、この各列の左側における行の入れ替えのキャンセルの説明図である。同図においてＬＵ分解され、作業領域Ｗに格納された結果が同じ大きさの作業領域Ｗ１にコピーされ、この作業領域Ｗ１上で、各列の左側で行われた行の入れ替えのキャンセルが行われる。列ブロックのＬＵ分解において行の入れ替えが行われた場合には、その入れ替えに関する情報はある１次元配列ＩＰ（ｎ）に格納されているものとし、その情報に基づいて左側における行の入れ替えのキャンセルが行われ、その結果作業領域Ｗ１の左下三角部分には再び０が格納されることになる。そしてこのキャンセル結果が圧縮モードの帯行列の格納領域、すなわち図４の対応する列に格納される。なお作業領域Ｗ上のＬＵ分解結果は、後述する行ブロックと行列の更新時に利用される。

この左側における行の入れ替えのキャンセルは、図４で説明した帯行列の格納形式と関連するものである。列ブロックのＬＵ分解に続いて対応する行ブロックと行列の更新が行われるが、左側で行われた行の入れ替え結果が残っていると、その後の行の入れ替えで、図４の下側の行数ｎｈ１の領域のさらに下側に対応する要素が０でなくなる可能性があり、メモリの格納領域を節約するために左側での行の入れ替えのキャンセルが行われる。なおここで更新対象として述べている行ブロックと行列については図１０で説明する。

その後、行ブロックの更新と、対応する行列の更新が行われるが、これらの更新に先立って更新によって壊れる可能性のある部分、すなわち図４の格納形式では本来０であるべき部分の要素が０でなくなる可能性のある領域のデータをあらかじめデータ退避用作業領域に退避しておき、その領域を０に設定した後に、行ブロックの更新と行列の更新が行われる。この部分は図６における退避必要領域である。すなわちこの三角部分は図４の配列形式において上の行数ｎｈ１の領域よりさらに上の領域に相当し、行の入れ替えによってこの部分が壊れる可能性がある場合には、この領域の退避が必要となる。

このように壊れる可能性のある領域の要素の値が退避された後に行ブロックと行列の更新が行われる。この更新においては、ＬＵ分解された列ブロックが圧縮モードの帯行列格納領域にコピーバックされた結果を用いて行われる。したがって行ブロックの更新、および行列の更新の対象となる行列は、圧縮モードで格納されている領域上での列のずれが補正された、行と列とが正常な位置関係を持つ配列として、例えば行ブロックの更新、および行列の更新を行うサブルーチンに渡されて、その後そのサブルーチンによる処理が行われる。この時、壊れる可能性があり、要素の値が退避された領域の要素に対しては、前述のようにその値が０に設定されて行の入れ替え、行ブロックの更新、および行列の更新が行われる。これらの更新において、前述の作業領域Ｗに残っている列ブロックのＬＵ分解結果が利用される。そして更新が終わった段階で、退避されていた壊れる可能性のある部分が元に戻される。

図１０は行ブロックの更新、および列ブロックと行ブロックに対応する行列の更新の並列処理の詳細説明図である。同図においてＵ１からＵ４までによって構成される行ブロックに対して、更新の前に列ブロックのＬＵ分解で行われた行の入れ替えに対応する入れ替えが行われる。この入れ替えも、例えばＵ１からＵ４、およびその下の行列のＣ１からＣ４を１次元目、すなわち行数によって分割し、４つのプロセッサにそれぞれ分担させて独立して並列処理することができる。

行ブロックの更新では、列ブロックの対角部分の下三角行列ＤＬの逆行列と行ブロックの分割部分のそれぞれＵｉとの行列積を並列に計算することによって、Ｕｉの更新処理が並列に行われる。そして残りの行列の更新では、分割された行列部分Ｃｉから、列ブロックの対角部分を除いたＬとＵｉとの積を減算することによって並列的にＣｉの更新処理が行われる。このように列ブロックのＬＵ分解から行列更新までの演算を並列処理することにより、メモリアクセス性能の低いスカラ計算機でも、メモリからロードしたデータに対する演算密度を高め、ＣＰＵの演算性能を引き出すことができる。

以上の処理によって帯行列のＬＵ分解が終了すると前進代入（前進消去）、および後退代入の処理によって連立方程式の解を求めるソルバー部分の処理が行われる。このソルバー部分の処理についても、従来においては前述のように帯行列が圧縮モードで格納され、ＬＵ分解の過程で行われた行の入れ替えに対して各列の左側の行の入れ替えがキャンセルされて格納されていたために、行列ベクトル積でベクトルの更新を行うことができなかった。このため本実施形態では、圧縮モードで格納されている帯行列を作業領域にコピーして処理を行う前に、列の左側でキャンセルされていた行の入れ替えを再び行った後に処理を実行するものとする。これによって行列ベクトル積を用いたベクトルの更新を行うことが可能となり、並列処理によって処理の効率を向上させることができる。

このソルバー部分の並列化処理について図１１、および図１２を用いて説明する。図１１は、帯行列から作業領域への列ブロックの格納（コピー）の説明図である。ＬＵ分解における処理と同様に、対角要素が最も左上にあるように複数の列から構成される列ブロックがある作業領域にコピーされる。前述と同様に左下の三角領域はもともと０の要素が格納されている領域であり、０クリアされる。そしてこの作業領域上で各列の左側に対する行の入れ替えが行われる。

この行の入れ替えは、ＬＵ分解の過程において行われた行の入れ替え情報、すなわち前述のＩＰ（ｎ）を用いて行われる。この入れ替え情報では、例えばｎ＝１行目に対して行われた入れ替えの相手側の行の情報が保存されている。例えばＩＰ（１）＝２２であれば１行目と２２行目との入れ替えが行われたことになり、また例えばＩＰ（１０）＝２であれば１０行目と２行目とが入れ替えられたことが入れ替え情報から判明し、この情報を使って各列の左側の行の入れ替えが行われる。

図１２は、このような列ブロックに対する行の入れ替えの後に、定数ベクトルｂの要素の入れ替えなどを行う前進消去における並列化処理の説明図である。同図においてまず行の入れ替えに対応して、入れ替え情報ＩＰ（ｎ）が保存されている、例えば１次元配列の情報を使って対応するベクトル、すなわち定数ベクトルｂ、および解ベクトルを格納する領域の入れ替えが行われ、列ブロックの対角部分の下三角行列に関する方程式が前進消去で解かれ、その後列ブロックのうちで対角部分の下の行列を用いて定数ベクトルの対応する部分の更新が行われる。この処理は並列的に行われる。

図１２において、４台のプロセッサを用いて並列処理を行うものとすれば、定数ベクトルｂを列ブロックの対角部分に対応するｂ０と、列ブロックの残りの部分の行数を均等に分割したｗ１からｗ４にそれぞれ対応する定数ベクトルの部分ｂ１からｂ４に分割した後に、対角部分の下三角行列ｗ１０に関する連立方程式ｗ１０×ｘ０＝ｂ０が解かれ、決定されたｘ０を用いてｂ０の更新が行われる。ここでｘ０は、解ｘのベクトルのうち、ベクトルｂ０に対応する部分である。

定数ベクトルの残りの部分ｂ１からｂ４に対しては、４台のプロセッサによって並列更新が行われる。すでに更新されたｂ０を用いて次式によってその更新処理が並列的に実行される。

ｂｉ＝ｂｉ−ｗｉ×ｂ０（ｉ＝１．．４）
以下本実施形態におけるＬＵ分解処理とソルバー処理とについてそれぞれフローチャートを用いてさらに詳細に説明する。図１３は、帯行列のＬＵ分解処理の詳細フローチャートである。同図において処理が開始されると、まずステップＳ１で行列の次数Ｎ、帯幅としての下バンド幅ｎｈ１、上バンド幅ｎｈ２、および帯行列が圧縮モードで格納された圧縮配列のデータが入力され、ステップＳ２でブロック幅ｎｂｌｋｓ、および繰返し数ｌｏｏｐが決定され、ステップＳ３でカウント数ｎｃｎｔの値が１とされる。ここで行列の次数Ｎが、例えば数十万であっても、ブロック幅ｎｂｌｋｓは例えば４０とされ、このブロック幅ｎｂｌｋｓを用いて繰返し数が次式によって決定される。

ｌｏｏｐ＝（Ｎ＋ｎｂｌｋｓ−１）／ｎｂｌｋｓ
すなわち繰返し数としては、単純に行列の次数Ｎをブロック幅ｎｂｌｋｓで割るのではなく、余りが出る場合を考えて、Ｎにｎｂｌｋｓ−１を加算してその結果をｎｂｌｋｓで割ったものが繰返し数とされる。

続いてステップＳ４で列ブロックの作業領域Ｗへのコピーが行われる。図３から図５で説明したように、まず最初に図３の最も左上の対角要素を含む最も左側の列から始め、ブロック幅ｎｂｌｋｓの数の列が作業領域Ｗにコピーされる。

このコピーされる列をカウント数ｎｃｎｔ、およびブロック数ｎｂｌｋｓを用いて一般的に表わすと、コピーされる列は（ｎｃｎｔ−１）×ｎｂｌｋｓ＋１列目からｎｃｎｔ×ｎｂｌｋｓ列目までであり、これらの列が作業領域Ｗ、すなわち（ｎｈ１＋ｎｂｌｋｓ）行、ｎｂｌｋｓ列の領域に行と列とが正常な位置関係となるようにコピーされる。

続いてステップＳ５で列ブロックのＬＵ分解と、その結果の作業領域Ｗへの格納が行われる。この列ブロックのＬＵ分解は前述のように特許文献１の方法を用いて並列処理として実行される。

列ブロックのＬＵ分解が終了すると、ステップＳ６でその結果が領域Ｗと同じ大きさを持つ作業領域Ｗ_１にコピーされ、ステップＳ７でその領域の上で各列の左側でピボット選択の結果として行われていた行の入れ替えがキャンセルされ、ステップＳ８でキャンセル結果が圧縮モードの帯行列の格納配列にコピーバックされる。

続いてステップＳ９で列ブロックのＬＵ分解の結果から、行の入れ替えを行う必要がある行の長さの最大値が計算され、ステップＳ１０でその最大値が圧縮モードの帯行列の配列領域を越える場合には、超えることによって壊れる可能性のある部分が作業領域に退避され、その部分が０に設定される。なおこの退避すべき領域は、圧縮モードの配列領域に格納された帯行列の行と列の関係を正常な位置関係に補正したものとして決定され、その領域の要素が０に設定される。

そしてステップＳ１１で圧縮モードの配列が実際に通常の位置関係、すなわち通常形式の行列に補正され、ステップＳ１２で行ブロックの並列更新において、行列の２次元目、すなわち列の数に応じて並列処理を行う各ＣＰＵの分担範囲が計算され、ステップＳ１３で行ブロックに対する行の入れ替えが列単位で並列に行われ、ステップＳ１４で作業領域Ｗに格納されている列ブロックのＬＵ分解の結果を用いて行ブロックの並列的更新が行われる。ステップＳ１５で列ブロックと行ブロックとを用いてこれらに対応する行列が更新され、ステップＳ１６で退避してあった部分が元に戻され、ステップＳ４で作業領域Ｗにコピーされた列ブロックに対応する処理を終了する。

そしてステップＳ１７でカウント数が繰返し数−１に達したか否かが判定され、達していない場合にはステップＳ１８でカウント数、すなわちｎｃｎｔの値がインクリメントされ、ステップＳ４以降の処理が繰り返される。

この処理では、まずカウント数ｎｃｎｔ＝１に対応して、例えば図３で最も左上に取られた列ブロックの対角部分に対応する行と列とを除いた小さくなった内側の行列部分に対してｎｃｎｔ＝２以上に対応する処理が繰り返され、ステップＳ１７でカウント数ｎｃｎｔが繰返し数−１に達したと判定されると、帯行列において最も右下の残った部分に対してステップＳ１９で１列ごとのＬＵ分解が行われて処理を終了する。

図１４は、ソルバー処理の詳細フローチャートである。同図において処理が開始されると、まずステップＳ２１からＳ２３において図１３のステップＳ１からＳ３と同様に各データの入力、ブロック幅と繰返し数の決定、およびカウント数の初期化が行われる。ステップＳ２４でステップＳ４におけると同様に、圧縮モードで格納されたＬＵ分解結果としての帯行列から列ブロックが作業領域Ｗにコピーされる。

そしてステップＳ２５で図１３の帯行列のＬＵ分解の過程で行われた行入れ替えの情報、例えば前述の１次元配列ＩＰ（ｎ）に保存されている入れ替え情報を用いて、作業領域Ｗ上で各列の左側の行の入れ替えが行われる。ステップＳ２６で定数ベクトルｂに対して対応する要素の入れ替えが行われ、また解ベクトルを格納する領域についても入れ替えが行われ、ステップＳ２７で列ブロックの対角部分の下三角行列に関する方程式が解かれ、ステップＳ２８で対角部分の下の行列が並列処理を行う各プロセッサによって分担され、行列ベクトル積を用いたベクトルの更新が行われる。ステップＳ２９でカウント数が繰返し数−１に達したか否かが判定され、まだ達していない場合にはステップＳ３０でカウント数がインクリメントされ、ステップＳ２４以降の処理が次の列ブロックに対応して実行される。

ステップＳ２９でカウント数が繰返し数−１に達したと判定されると、ステップＳ３１で最後の対角ブロックに関する連立方程式の解が求められる。この処理では、行の入れ替え情報を用いてベクトルの要素を入れ替えた後で、残りのベクトル部分を更新する処理が最後の要素まで繰り返される。

その後ステップＳ３２で後退代入の計算が行われる。後退代入の処理においては、行の入れ替えの必要はなく、ブロック幅ごとに最後、すなわち下の方から順番に処理が行われる。この処理では対角部分の上三角行列に関する連立方程式が解かれ、解のベクトルを利用してそれより上の部分のベクトルが行列ベクトル積で更新される処理が行われ、後退代入の計算が終了すると、ソルバー処理を終了する。

以上において、本実施形態における連立１次方程式の解法プログラムについて詳細に説明したが、このプログラムを用いることにより帯係数行列を持つ連立１次方程式の計算が高速化される。１つの例として従来のベクトル計算機向けのコードと比較した場合、１ＣＰＵのスカラ計算機でも１４倍高速となることが判明した。共用メモリ型スカラ計算機を使用することにより、さらなる高速化が期待できる。

（付記１） 帯係数行列を持つ連立方程式の解法プログラムであって、
前記帯行列内の複数の列から構成される列ブロックを並列的にＬＵ分解し、作業領域に格納する手順と、
該ＬＵ分解結果に対して各列の左側でピポット選択の結果として行われていた行の入れ替えをキャンセルし、該キャンセル結果のデータを圧縮形式の帯行列格納配列領域にコピーバックする手順と、
前記列ブロックのＬＵ分解結果に対応して、入れ替えの必要な行の長さの最大値と前記圧縮形式の帯行列格納配列領域の大きさとに対応して、帯行列の更新によって壊れる可能性のある部分のデータを退避する手順と、
前記作業領域に格納されている前記列ブロックのＬＵ分解結果を用いて、前記帯行列を並列的に更新する手順と、
前記退避部分のデータを該帯行列更新結果上に戻す手順とを含む行列演算を計算機に実行させることを特徴とする帯係数行列を持つ連立方程式の解法プログラム。

（付記２） 前記行列演算が、前記帯行列のＬＵ分解のための演算であることを特徴とする付記１記載の帯係数行列を持つ連立方程式の解法プログラム。
（付記３） 前記帯行列のＬＵ分解のための演算において、
前記帯行列の最も左上部の前記列ブロックに対応する前記行列演算の終了後に、該列ブロック最上部の対角部分に含まれる行と列とを前記帯行列から除外し、該除外後の行列の左上部の列ブロックに対応する前記行列演算を繰返し、最後に残った部分のＬＵ分解を行うことを特徴とする付記２記載の帯係数行列を持つ連立方程式の解法プログラム。

（付記４） 前記最後に残った部分のＬＵ分解の終了後に、前記帯行列の並列的更新において行われた行入れ替えの情報を利用して、複数の列から構成される列ブロックの各列の左側で、ピボット選択に対応する行の入れ替えを行うとともに、該左側の行の入れ替えに対応して連立方程式における定数ベクトルの要素の内で該列ブロックに対応する要素の入れ替えを行う手順と、
該列ブロックの対角部分の下三角行列に対応する連立方程式を解く手順と、
該列ブロックの対角部分の下の行列を用いて行列ベクトル積の計算によってベクトルの更新を行う手順とを含む並列的行列演算をさらに計算機に実行させることを特徴とする付記３記載の帯係数行列を持つ連立方程式の解法プログラム。

（付記５） 前記並列的行列演算が、前記帯行列のＬＵ分解結果を用いた前進代入処理の演算であり、
前記帯行列の最も左上部の前記列ブロックに対応する前記並列的行列演算の終了後に、該列ブロック最上部の対角部分に含まれる行と列とを前記帯行列から除外して、該除外後の行列の左上部の列ブロックに対応する該並列的行列演算を繰返し、最後に残った対角ブロックの下三角行列に対する連立方程式の解を求めることを特徴とする付記４記載の帯係数行列を持つ連立方程式の解法プログラム。

（付記６） 前記帯行列の並列的更新手順において、
前記キャンセル結果のデータがコピーバックされ、圧縮形式の配列領域に格納された帯行列に対して、各列の要素の行の位置を補正して補正後の行列の更新を行うことを特徴とする付記１記載の帯係数行列を持つ連立方程式の解法プログラム。

（付記７） 前記帯行列の並列的更新手順において、
前記補正後の帯行列と前記列ブロックのＬＵ分解結果とを用いて、該列ブロック内の対角行列に対応する行ブロックを並列的に更新する手順と、
該行ブロックの更新結果を用いて該列ブロックと行ブロックとに対応する行列を並列的に更新する手順とを備えることを特徴とする付記６記載の帯係数行列を持つ連立方程式の解法プログラム。

（付記８） 前記連立方程式の解法プログラムが共用メモリ型スカラ計算機によって実行されることを特徴とする付記１記載の帯係数行列を持つ連立方程式の解法プログラム。
（付記９） 帯係数行列を持つ連立方程式の解法プログラムであって、
該帯係数行列のＬＵ分解演算の終了後に、該帯係数行列の更新において行われた行入れ替えの情報を利用して、複数の列から構成される列ブロックの各列の左側で、ピボット選択に対応する行の入れ替えを行うとともに、該左側での行の入れ替えに対応して連立方程式における定数ベクトルの要素の内で該列ブロックに対応する要素の入れ替えを行う手順と、
該列ブロックの対角部分の下三角行列に対応する連立方程式を解く手順と、
該列ブロックの対角部分の下の行列を用いて、行列ベクトル積の計算によってベクトルの更新を行う手順とを含む並列的行列演算を計算機に実行させることを特徴とする帯係数行列を持つ連立方程式の解法プログラム。

（付記１０） 前記連立方程式の解法プログラムが共用メモリ型スカラ計算機によって実行されることを特徴とする付記９記載の帯係数行列を持つ連立方程式の解法プログラム。

本発明の連立１次方程式解法プログラムの原理的な機能ブロック図である。 本発明のプログラムを並列処理として実行する共用メモリ型スカラ並列計算機のハードウェア構成例を示す図である。 帯行列の形式を説明する図である。 帯行列の圧縮モード格納形式の説明図である。 圧縮モードで格納された帯行列から列ブロックの作業領域への格納の説明図である。 列ブロックの帯行列内の位置を説明する図である。 列ブロックの下三角部分への０設定の説明図である。 列ブロックのＬＵ分解における並列処理の説明図である。 列ブロックのＬＵ分解結果における各列の左側の行の入れ替えのキャンセルの説明図である。 行ブロックと行列の並列更新処理の説明図である。 ソルバー部分の並列解法における列ブロックの作業領域への格納の説明図である。 列ブロックに対応する定数ベクトルの並列更新処理の説明図である。 帯行列のＬＵ分解処理の詳細フローチャートである。 ソルバー処理の詳細フローチャートである。

符号の説明

１０ プロセッサ
１１ メモリモジュール
１２ 相互結合網
１３ ２次キャッシュメモリ

Claims (5)

1. 帯係数行列を持つ連立方程式の解法プログラムであって、
   前記帯行列内の複数の列から構成される列ブロックを並列的にＬＵ分解し、作業領域に格納する手順と、
   該ＬＵ分解結果に対して各列の左側でピポット選択の結果として行われていた行の入れ替えをキャンセルし、該キャンセル結果のデータを圧縮形式の帯行列格納配列領域にコピーバックする手順と、
   前記列ブロックのＬＵ分解結果に対応して、入れ替えの必要な行の長さの最大値と前記圧縮形式の帯行列格納配列領域の大きさとに対応して、帯行列の更新によって壊れる可能性のある部分のデータを退避する手順と、
   前記作業領域に格納されている前記列ブロックのＬＵ分解結果を用いて、前記帯行列を並列的に更新する手順と、
   前記退避部分のデータを該帯行列更新結果上に戻す手順とを含む行列演算を計算機に実行させることを特徴とする帯係数行列を持つ連立方程式の解法プログラム。
2. 前記行列演算が、前記帯行列のＬＵ分解のための演算であることを特徴とする請求項１記載の帯係数行列を持つ連立方程式の解法プログラム。
3. 前記帯行列のＬＵ分解のための演算において、
   前記帯行列の最も左上部の前記列ブロックに対応する前記行列演算の終了後に、該列ブロック最上部の対角部分に含まれる行と列とを前記帯行列から除外し、該除外後の行列の左上部の列ブロックに対応する前記行列演算を繰返し、最後に残った部分のＬＵ分解を行うことを特徴とする請求項２記載の帯係数行列を持つ連立方程式の解法プログラム。
4. 前記最後に残った部分のＬＵ分解の終了後に、前記帯行列の並列的更新において行われた行入れ替えの情報を利用して、複数の列から構成される列ブロックの各列の左側で、ピボット選択に対応する行の入れ替えを行うとともに、該左側の行の入れ替えに対応して連立方程式における定数ベクトルの要素の内で該列ブロックに対応する要素の入れ替えを行う手順と、
   該列ブロックの対角部分の下三角行列に対応する連立方程式を解く手順と、
   該列ブロックの対角部分の下の行列を用いて行列ベクトル積の計算によってベクトルの更新を行う手順とを含む並列的行列演算をさらに計算機に実行させることを特徴とする請求項３記載の帯係数行列を持つ連立方程式の解法プログラム。
5. 帯係数行列を持つ連立方程式の解法プログラムであって、
   該帯係数行列のＬＵ分解演算の終了後に、該帯係数行列の更新において行われた行入れ替えの情報を利用して、複数の列から構成される列ブロックの各列の左側で、ピボット選択に対応する行の入れ替えを行うとともに、該左側での行の入れ替えに対応して連立方程式における定数ベクトルの要素の内で該列ブロックに対応する要素の入れ替えを行う手順と、
   該列ブロックの対角部分の下三角行列に対応する連立方程式を解く手順と、
   該列ブロックの対角部分の下の行列を用いて、行列ベクトル積の計算によってベクトルの更新を行う手順とを含む並列的行列演算を計算機に実行させることを特徴とする帯係数行列を持つ連立方程式の解法プログラム。

Images