(57)【要約】【課題】共有メモリ型スカラ計算機において、高速に固有値問題の演算が行える演算の並列処理方法を提供する。【解決手段】固有値問題を解く方法は、３段階に分かれ、行列を３重対角化するステップと、３重対角行列から固有値と固有ベクトルを求めるステップと、３重対角行列から求まった固有ベクトルを変換して、元の行列の固有ベクトルを求めるステップとからなる。特に、３重対角化を行うステップと、元の行列の固有ベクトルを求めるステップの演算負荷が大きいので、これらのステップを並列化して処理を行い、高速に固有値問題を解く。【選択図】　　　図２
2004005528.tif 000002
【特許請求の範囲】
【請求項１】
共有メモリ型スカラ並列計算機用固有値問題の並列処理方法を実現させるプログラムであって、
演算すべき実対称行列あるいはエルミート行列をブロック化し、ブロック化した部分をメモリの作業域にコピーして、ブロック間の積を利用して３重対角化するステップと、
３重対角化された行列から固有値と固有ベクトルを求めるステップと、
該３重対角化された行列から求められた固有ベクトルを、所定のブロック幅の行列演算の並列演算となるように変形されたハウスホルダー変換によって変換し、元の行列の固有ベクトルを求めるステップと、
を備えることを特徴とする並列処理方法を共有メモリ型スカラ並列計算機に実現させるプログラム。
【請求項２】
前記３重対角化ステップにおいて、ブロック化された部分の更新を再帰的プログラムで行うことを特徴とする請求項１に記載のプログラム。
【請求項３】
前記３重対角化ステップにおいて、ブロック化された部分は、キャッシュメモリのページを跨らないように小さなブロックに分割して格納し、各プロセッサで並列に演算を行うことによって処理されることを特徴とする請求項１に記載のプログラム。
【請求項４】
前記固有ベクトルを求めるステップにおいて、前記ハウスホルダー変換のための行列は、各プロセッサで処理できるような小さな係数行列となる上三角行列を各プロセッサで同時に生成されることを特徴とする請求項１に記載のプログラム。
【請求項５】
前記固有ベクトルを求めるステップにおいて、前記固有ベクトルは、格納されている２次元配列の２次元目をプロセッサの数に対応して均等に分割し、おのおのの分割された領域を各プロセッサに割り当てることによって、演算されることを特徴とする請求項１に記載のプログラム。
【発明の詳細な説明】
【０００１】
【発明の属する技術分野】
本発明は、共有メモリ型スカラ並列計算機における行列の演算に関する。
【０００２】
【従来の技術】
ベクトル並列計算機用実対称行列（行列要素が実数で、転置しても行列が不変に保たれる行列）及びエルミート行列（行列要素が複素数で、共役・転置をとっても不変に保たれる行列）の固有値問題（Ａを行列、λを定数、Ｉを単位行列とした場合、ｄｅｔ｜Ａ−λＩ｜＝０となるλとその固有ベクトルを求めること）は、３重対角化（対角成分とその両隣の成分のみを持つ行列にすること）をまず行う。この３重対角行列の固有値問題をマルチセクション法によって解き、固有値を求めてから逆反復法で固有ベクトルを求める。その後、固有ベクトルにハウスホルダー変換を行い、元の固有値問題の固有ベクトルを求めていた。
【０００３】
【発明が解決しようとする課題】
ベクトル計算機においては、メモリアクセス系の動作速度が速いことを前提として、固有値問題の演算を行っていたが、共有メモリ型スカラ計算機の場合、演算する行列が大きくなるほどに、共有メモリにアクセスする回数が多くなり、アクセス速度の遅い共有メモリへのアクセスによって、計算機の性能が大きく損なわれてしまうという問題がある。そこで、共有メモリ型スカラ計算機の各プロセッサに設けられる、アクセス速度の速いキャッシュメモリを有効に使い、行列計算を行う必要がある。つまり、行列の各行あるいは列毎に演算していると、共有メモリへのアクセスが多くなってしまうので、行列をブロック化して、各プロセッサにキャッシュメモリに格納されたデータを最大限処理した後、共有メモリにアクセスするようにして、共有メモリへのアクセス数を減少する、各プロセッサに局所化したアルゴリズムが必要となる。
【０００４】
すなわち、共有メモリ型の並列計算機では、ベクトル計算機のような強力なメモリアクセス能力がないため、メモリアクセスに対して演算量を増やすアルゴリズム上の工夫が必要になる。
【０００５】
本発明の課題は、共有メモリ型スカラ計算機において、高速に固有値問題の演算が行える演算の並列処理方法を提供することである。
【０００６】
【課題を解決するための手段】
本発明の並列処理方法は、共有メモリ型スカラ並列計算機用固有値問題の並列処理方法を実現させるプログラムであって、演算すべき実対称行列あるいはエルミート行列をブロック化し、ブロック化した部分をメモリの作業域にコピーして、ブロック間の積を利用して３重対角化するステップと、３重対角化された行列から固有値と固有ベクトルを求めるステップと、該３重対角化された行列から求められた固有ベクトルを、所定のブロック幅の行列演算の並列演算となるように変形されたハウスホルダー変換によって変換し、元の行列の固有ベクトルを求めるステップとを備えることを特徴とする。
【０００７】
本発明によれば、固有値問題を共有メモリ型スカラ計算機の各プロセッサにできるだけ局所化した演算で解くことができるので、各プロセッサが、共有メモリに頻繁にアクセスするために生じる遅延時間が最小限に抑えられ、並列計算の威力を有効に発揮した演算が行える。
【０００８】
【発明の実施の形態】
本発明の実施形態においては、上記固有値問題の演算における３重対角化部分のアルゴリズムでは、ブロック化されたアルゴリズムを採用し、ブロック部分を求める部分のアルゴリズムを再帰的に行うようにし、更新部分での演算密度を高める。また、メモリのアクセスがキャッシュメモリのページ間をとびとびにアクセスすることを避けるように行列ベクトル積部分を対称性を利用して連続アクセスが使える方式にする。キャッシュメモリのページを跨ってデータを読み出す場合、一度に読み込めず、キャッシュメモリに２度アクセスする必要が出てくるので、計算機のパフォーマンスを落とすことにある。従って、上記のように、ページをデータがまたがないようにする。
【０００９】
また、３重対角化された行列の固有ベクトルにハウスホルダー変換を行い、元の行列の固有ベクトルを求める部分は、ハウスホルダー変換を８０回程度ごとに束ねて、３つの行列要素からなる演算を行うことで演算密度を高める。
【００１０】
なお、本発明の実施形態においては、３重対角化された行列から固有値を求める方法及び３重対角化された行列の固有ベクトルを求める方法は、従来の方法によるとする。
【００１１】
図１は、本発明の実施形態が前提とする共有メモリ型スカラ計算機のハードウェア構成を示した図である。
プロセッサ１０−１〜１０−ｎは、１次キャッシュメモリを持っており、この１次キャッシュメモリはプロセッサの中に組み込まれていることもある。また、各プロセッサ１０−１〜１０−ｎには、２次キャッシュメモリ１３−１〜１３−ｎが設けられ、２次キャッシュメモリ１３−１〜１３−ｎが相互結合網１２に結合されている。また、相互結合網１２には、共有メモリであるメモリモジュール１１−１〜１１−ｎが設けられ、プロセッサ１０−１〜１０−ｎは、演算に必要なデータをここから読み出し、相互結合網１２を介して、２次キャッシュメモリ１３−１〜１３−ｎあるいは、１次キャッシュメモリに記憶させて、演算を行う。
【００１２】
この場合、メモリモジュール１１−１〜１１−ｎから２次キャッシュメモリ１３−１〜１３−ｎあるいは１次キャッシュメモリにデータを読み込んだり、２次キャッシュメモリ１３−１〜１３−ｎ、あるいは、１次キャッシュメモリから演算後のデータをメモリモジュール１１−１〜１１−ｎに書き込むのはプロセッサ１０−１〜１０−ｎの演算速度に比べて非常に遅い。従って、このような書き込み、読み出しが頻繁に発生すると計算機全体の性能を劣化させてしまう。
【００１３】
従って、計算機全体の性能を高く維持するためには、メモリモジュール１１−１〜１１−ｎへのアクセスをできるだけ減らし、２次キャッシュメモリ１３−１〜１３−ｎ、１次キャッシュメモリ、及びプロセッサ１０−１〜１０−ｎからなるローカルな系においてできるだけ多くの演算をしてしまうようなアルゴリズムが必要となる。
＜固有値、固有ベクトルを求める方式＞
１．３重対角化部分
１）３重対角化
ａ）ブロック化された３重対角化の数学的アルゴリズム
行列をブロック幅ごとに３重対角化する。行列をブロックに分割してブロック単位で下記のアルゴリズムで３重対角化を行う。
【００１４】
図２〜図１１は、本発明の実施形態のアルゴリズムを説明する図である。
図２は、ｍ番目のブロックを処理しているときの図である。ブロックとは図２で点線で示した列及び行方向の矩形である。
【００１５】
最後のブロックは、そのブロック幅−２をブロック幅とする、右下に２×２行列を残した左詰めのブロックに関してアルゴリズムを適用して全ての処理を終わる。
【００１６】
ｄｏ　ｉ＝１，ｂｌｋｓ
ｓｔｅｐ１：Ａｎのｎ＋１行ベクトルよりハウスホルダーベクトルｕを作る。
ｓｔｅｐ２：ｖｉ＝Ａｎ＋ｉｕ、ｗｉ＝ｖｉ−ｕ（ｕｔｖ）／２を計算する。
ｓｔｅｐ３：Ｕｉ＝（Ｕｉ−１、ｕｉ）、Ｗｉ＝（Ｗｉ−１、ｗｉ）と更新する。（ここで、（Ｕｉ−１、ｕｉ）は、行列Ｕｉ−１という１列を加えて行列Ｕｉとすることにより、行列を１列づつ拡大する。）
ｓｔｅｐ４：　ｉｆ（ｉ＜ｂｌｋｓ）　ｔｈｅｎ
Ａｎのｎ＋ｉ＋１列目を更新する。
【００１７】
2004005528.tif 000003
・ブロック化されたハウスホルダー変換による３重対角化
・・ハウスホルダー変換の説明
ｖ＝（ｖ１，ｖ２，．．．．．．，ｖｎ）
｜ｖ｜２＝ｖ＊ｖ＝ｈ２
Ｕｖ＝（ｈ，０，．．．．，０）ｔとすると、Ｕｖ＝ｖ−（ｖｉ−ｈ，ｖ２，．．．．．，ｖｎ）の関係がある。
Ｕ＝（１−ｕｕｔ／｜ｕ｜２）＝（１−αｕｕｔ）、ここでｕ＝（ｖ１−ｈ，ｖ２，．．．．，ｖｎ）
以下、αを省略して計算する。
2004005528.tif 000004
これを繰り返し利用すれば、
Ａｎ＋ｋ＝Ａｎ−ＵｋＷｋｔ−ＷｋＵｋｔ　　　　　．．．．．．．．．．．（＊＊）
ｋステップ目の計算は（＊）及び（＊＊）からＶｎを次のように計算できる。
ｖｋ＝Ａｎｕｋ−Ｕｋ−１Ｗｋ−１ｔｕｋ−Ｗｋ−１Ｕｋ−１ｔｕｋ　　　　．．．．．．．．（＊＊＊）
ｗｋ＝ｖｋ−ｕｋｕｋｔｖｋ／２
Ｕｋ＝（Ｕｋ−１，　ｕｋ），　Ｗｋ＝（Ｗｋ−１，　ｗｋ）
Ａｎ＋ｋ＝Ａｎ−ＵｋＷｋｔ−ＷｋＵｋｔ
ｂ）ハウスホルダー変換を構成する情報の格納
固有ベクトルを計算するときに３重対角化で使ったハウスホルダー変換が必要になる。このためＵｎおよびαをハウスホルダー変換を構成するベクトルの位置に格納する。αは対応する対角要素の位置に格納する。
ｃ）Ｕｉを効率的に求める方法
ブロック部分の３重対角化を行うためにハウスホルダー変換のための次のベクトルを更新する必要がある。これらの計算をできるだけローカルに行うために、ブロック幅部分を作業域にコピーして３重対角化を行い、元の領域に格納する。次の列ベクトルの更新をそのたび行わず、行列積の形にして演算密度を上げて計算を行う。このため、ブロック部分の３重対角化を再帰的なプログラミングで行う。
【００１８】
ｒｅｃｕｒｓｉｖｅ　ｓｕｂｒｏｕｔｉｎｅ　ｔｒｉｄ（幅、ブロック領域ポインタ）
ｉｆ（幅＜１０）ｔｈｅｎ
ｃ　当該幅のブロック部分の３重対角化を行う。
【００１９】
ハウスホルダー変換で必要なベクトルｕ、及び行列ベクトル積からｖｉおよびｗｉを作成する。
Ｕ、Ｗにｕｉ、ｗｉを結合する。
【００２０】
ｅｌｓｅ
ｃ　ブロック幅を半分に分割
ｃ　前半ブロックに関して３重対角化を行う。
【００２１】
ｃａｌｌ　ｔｒｉｄ（前半の幅、前半の領域）
ｃ　ブロックを分割して、分割線から後半部分を更新
Ｂ＝Ｂ−ＵＷｔ−ＷＵｔで更新する。
ｃ　その後、後半部分を３重対角化する。
【００２２】
ｃａｌｌ　ｔｒｉｄ（後半の幅、後半の領域）
ｒｅｔｕｒｎ
ｅｎｄ
図３に示されるように、作業領域Ｕにブロック部分をコピーして、ブロック部分を再帰的プログラムで３重対角化する。再帰的プログラムのため、図３の前半部分は前半の更新処理のために当該再帰的プログラムが呼び出されたときに３重対角化される。後半部分は、前半部分によって更新され、その後、３重対角化される。
【００２３】
図４に示されるように、再帰的プログラムが深さ２まで呼ばれると、最初の前半処理で図４（ａ）の斜線部分をＢとして更新し、次に、図４（ｃ）の斜線部分、最後に図４（ｅ）が更新される。更新時の並列化は、更新部分のブロック行列を行ベクトル方向に均等に分割して、各部分の更新を複数のプロセッサで並列に行う。
図４（ａ）の演算後、図４（ｂ）の演算が行われ、図４（ｃ）の演算後、図４（ｄ）の演算が行われ、図４（ｅ）の演算後、図４（ｆ）の演算が行われる。
【００２４】
図５に示されるように、Ｕの斜線部分を更新するとき、ｕの横線及びＷの縦線部分を参照する。このことで、演算密度を上げることができる。すなわち、
Ａｎ＋ｋ＝Ａｎ−ＵｋＷｋｔ−ＷｋＵｋｔ　　　．．．．．．．（＊＊）
を利用して、Ｖｎを次のように計算できる。このときＵ、Ｗの参照パターンは
ｖｋ＝Ａｎｕｋ−Ｕｋ−１Ｗｋ−１ｔｕｋ−Ｗｋ−１Ｕｋ−１ｔｕｋ　　　　．．．．．（＊＊＊）
図５（ａ）と（ｂ）、（ｃ）と（ｄ）、（ｅ）と（ｆ）のＵの更新の後、更新した部分の３重対角化のためにｖｋを計算するときには、Ｕ、Ｗを参照して、行列ベクトル積で計算する。参照だけであり、Ｕの更新と参照が共通する部分もあり、効率的に参照できる。Ａｎを毎回更新せずに、このようにＵ、Ｗを使って必要な部分を更新し、（＊＊）を使って、全体の更新の演算密度を上げて性能向上を行う。（＊＊＊）は余分な演算ではあるが、ブロック幅を小さくしておく限り無視できる。
【００２５】
例えば、４個のＣＰＵで並列処理を行う場合、（＊＊＊）のＷｋ−１ｔｕｋ、Ｕｋ−１ｔｕｋの計算では、縦線での分割で並列計算を行い、その結果との積に関しては、破線の方向に分割して並列計算する。
ｖｉ＝Ａｎｕｉの計算の並列化
図６に示されるように、Ａｎの対称性、つまり、Ａｎ＝Ａｎｔを利用して各プロセッサで２次元目を分割して、ｖｉの部分をＡｎ（＊、ｎｓ：ｎｅ）ｔｕｉを各プロセッサで計算する。
２）共有メモリ型スカラ並列計算機での並列化
ａ）Ｕ、Ｗ用の格納領域は、共有メモリに確保する。また、３重対角化を行うブロック領域は別に確保した作業領域にコピーして、その領域で３重対角化を行う。
【００２６】
上記の再帰的プログラムの中で、並列化する部分は以下の部分
（１）ハウスホルダー変換を構築するために必要なｕｉを求めるためにｓｔｅｐ４の
Ａｎ（＊，　ｎ＋ｉ＋１）＝Ａｎ（＊，　ｎ＋ｉ＋１）−ＵｉＷｉ（ｎ＋ｉ＋１，＊）ｔ−ＷｉＵｉ（ｎ＋ｉ＋１，＊）ｔ
を使って、必要なベクトルを求める部分
（２）ｓｔｅｐ２でｖｉを計算する部分。この部分は（＊＊）の以下の式にｕｉを作用させて計算する。
Ａｎ＋ｋ＝Ａｎ−ＵｋＷｋｔ−ＷｋＵｋｔ
この計算ではＡｎとｕｉとの積及びＵｋＷｋｔ−ＷｋＵｋｔとｕｉとの積を並列化する。
【００２７】
ブロック部分は、作業域にコピーされておりＡｎの必要部分は更新されないようにする。また、Ａｎの対称性を利用して列ベクトル方向に長い行列に分割して並列計算する。
（３）再帰的プログラムで、ブロック領域の更新は、
Ａｎ＋ｋ＝Ａｎ−ＵｋＷｋｔ−ＷｋＵｋｔ
を利用して更新する。そのことで（１）の計算を減らす。
３）ｓｔｅｐ５の更新部分
更新は対称性を利用して対角成分より下半分を計算する。並列化は負荷がバランスするようにＣＰＵ数を＃ＣＰＵとしたとき、更新すべき部分行列の格納されているサブアレイの２次元目を２×＃ＣＰＵ個に均等に分割して、１から２×＃ＣＰＵと番号を振る。１から＃ＣＰＵのＣＰＵに関してｉ番目のプロセッサでｉ番目の部分と２×＃ＣＰＵ＋１−ｉ番目の部分の更新を並列に行う。
【００２８】
このあと上半分にコピーを行う。これも同様に分割して負荷をバランスさせる。このとき、対角ブロック部分以外は、キャッシュメモリのページをまたがらないようにある程度小さなブロックに分割してコピーを行う。
Ａｎ＋ｋ＝Ａｎ−ＵｋＷｋｔ−ＷｋＵｋｔで更新する部分は下三角行列部分で、図７のように列ブロックにＣＰＵ数×２に分割して、両端からペアを作り、これを各ＣＰＵで更新する。図７は、４個のＣＰＵが設けられている場合を記述している。
【００２９】
下三角部分が更新できたら、上三角部分に１〜８を同じペアで転置してｕ１からｕ８にコピーする。
このとき、小さな正方ブロックに内部を分割してキャッシュを利用して転置する。更新と同様にＣＰＵで並列処理する。
キャッシュ上での転置で性能向上する部分の説明
図８に示されるように、１から番号順に転置変換していく。各正方形領域に関して、１の場合は下三角部分をメモリの連続領域にコピーして、それを行方向にアクセスして（転置して）１の上三角部分に格納する。２〜８は正方形にコピーして、行方向に参照して列方向（上下方向）に格納領域に格納する。
２固有ベクトルの計算
ａ）基本的なアルゴリズム
ベクトルｕｎを格納して、後で（１−２＊ｕｕｔ／（ｕｔｕ））を構成して、ベクトルにかける。
【００３０】
３重対角化を行うと元の固有値問題は以下のようになる。
Ｑｎ−２．．．．．．Ｑ２Ｑ１ＡＱ１ｔＱ２ｔ．．．．．．．．Ｑｎ−２ｔＱｎ−２．．．．．．Ｑ２Ｑ１ｘ＝λＱｎ−２．．．．．．Ｑ２Ｑ１ｘ
３重対角化した固有値問題を解いて求めた固有ベクトルｙからｘ＝Ｑ１ｔＱ２ｔ．．．．．Ｑｎ−３ｔＱｎ−２ｔｙを計算することで変換することができる。
ｂ）本発明の実施形態のブロックアルゴリズム及び固有ベクトルの変換の並列計算
多数のまたは全ての固有ベクトルを求める場合は各ＣＰＵに均等に３重対角行列の固有ベクトルを対応付けて各ＣＰＵで並列に上記の変換を施す。このとき変換行列を８０本程度づつまとめて変換を行う。
【００３１】
各々の変換行列Ｑｉｔは１＋αｉｕｉｕｉｔと書ける。これらの積を展開すると以下の形になる。
ｂｉ，ｊは最左端及び最右端のｕｉｕｊを除いたスカラ係数をまとめたものである。
【００３２】
【数１】
2004005528.tif 000005
【００３３】
の形になる。ｂｉ，ｊは上三角行列になる。これを求めて１＋ＵＢＵｔの形になる。このような形に変形することで、計算密度を向上させることで高速化を図る。図９は、行列Ｂの模式図である。
【００３４】
上記の方法は３ステップであるが、メモリアクセスの関係から行列はＵ及びＢの２つとなり、Ｂは比較的小さな領域にすることができるため効率がよい。
ｍ−１番目までｂｉ，ｊができた後に、（１＋αｍＵｍＵｍｔ）をかけると、
【００３５】
【数２】
2004005528.tif 000006
【００３６】
最後の項の和をｉ、ｊを入れ替えて
【００３７】
【数３】
2004005528.tif 000007
【００３８】
と変更することができる。
最内の括弧の中がｂｍ，ｊ（ｊ＝ｍ＋１、・・・、ｎ＋ｋ）と考えられる。ｂｍ，ｍはαｍである。
ブロック幅の正方な作業配列Ｗ２を用意して、はじめにα１ＵｉＵｊｔをｗ２（ｉ、ｊ）の上三角部分に格納する。対角要素にはαｉを格納する。
【００３９】
次に、以上の議論は最右下の２×２上三角行列から順番に１行ずつ上に付け加える要領で計算することができる。
行の最も右の要素から一つずつ計算していけば、Ｂは上三角行列であるため更新された部分の参照が無いために同じ領域で計算ができる。このように、極めて小さな領域のみを使って３つの行列積の真ん中に位置する係数行列を計算することができる。
【００４０】
図１０は、上記固有値を求める方法を模式的に示した図である。
ブロック幅をｎｂｓとする。
まず、Ｂの上半分に内積αｊｕｉ・ｕｊを計算して格納する。
【００４１】
対角要素にはαｉを格納する。
次に、以下のように計算する。
ｄｏ　ｉ１＝ｎｂｓ−２，１，−１
ｄｏ　ｉ２＝ｎｂｓ，ｉ１＋１，−１
ｓｕｍ＝ｗ２（ｉ１，ｉ２）
ｄｏ　ｉ３＝ｉ２−１，　ｉ１＋１，　−１
ｓｕｍ＝ｓｕｍ＋ｗ２（ｉ１，ｉ３）＊ｗ２（ｉ３，　ｉ２）
ｅｎｄｄｏ
ｗ２（ｉ１，ｉ２）＝ｓｕｍ
ｅｎｄｄｏ
ｅｎｄｄｏ
ｄｏ　ｉ２＝ｎｂｓ，１，−１
ｄｏ　ｉ１＝ｉ２−１，１，−１
ｗ２（ｉ１，ｉ２）＝ｗ２（ｉ１，ｉ２）＊ｗ２（ｉ２，ｉ２）
ｅｎｄｄｏ
ｅｎｄｄｏ
図１１は、上記で求められた固有ベクトルを元の行列の固有ベクトルに変換する様子を模式的に示した図である。
【００４２】
配列Ａに格納されているハウスホルダーベクトルで固有ベクトルを変換する。ブロック化されていて、図１１の斜線の部分をＥＶの斜線部分にかけて結果をＷに格納する。また、Ａブロック行列からＷ２を構成する。Ｗ２とＷをかけたものにＡのブロック部分をかけたものからＥＶの斜線部分を更新する。
３エルミート行列の固有値・固有ベクトル
エルミート行列の固有値・固有ベクトルを求めるアルゴリズムは実対称行列の３重対角化での転置を転置＋複素共役に置き換える（ｔ→Ｈ）。また、ハウスホルダーベクトルの作り方は最初の要素のスカラ倍になるようにベクトルの大きさを加える。
【００４３】
できた３重対角行列はエルミート行列であり、これを絶対値１の対角行列でスケーリングする。
対角行列は以下のように作る。
ｄｉ＝１．０、ｄｉ＋１＝ｈｉ＋１，ｉ／｜ｈｉ＋１｜＊ｄｉ
図１２〜図１８は、本発明の実施形態に従ったルーチンの疑似コードである。
【００４４】
図１２は、実対称行列を３重対角化するサブルーチンである。
配列ａに実対称行列の下三角部分を格納する。ｄａｉｇ、ｓｄｉａｇに３重対角行列、副対角部分を格納する。出力としてａの下三角部分に変換で必要な情報を格納する。
【００４５】
Ｕは３重化を行うブロックを格納する。ＶはＷを格納する領域である。
ｎｂは、ブロックの数であり、ｎｂａｓｅはブロックの始まり位置を示す。
サブルーチンｃｏｐｙを実行した後、３重対角化をするブロックをｕ（ｎｂａｓｅ＋１：ｎ、１：ｉｂｌｋ）に格納し、ｂｌｋｔｒｉｄルーチンを呼び出してＬＵ分解を行う。そして、処理の終わったｕ（ｎｂａｓｅ＋１：ｎ、１：ｉｂｌｋ）を元の行列ａに書き戻す。これらの下の処理は、最後の残りのブロックの３重対角化をサブルーチンｂｌｋｔｒｉｄを用いて行うものである。
【００４６】
図１３は、３重対角化を行うサブルーチンの疑似コードである。
このサブルーチンは、ブロック行列を３重対角化するルーチンで、再帰的に呼び出される。ｎｂａｓｅはブロックの位置を示すオフセット、ｉｓｔａｒｔは再帰呼び出しで対象となる縮小されたサブブロックのブロック内でのオフセットで最初１で再帰的に呼び出されるとき対象サブブロックの位置を示す。ｎｗｉｄｔｈはサブブロックの大きさである。
【００４７】
ｎｗｉｄｔｈが１０より小さい場合には、サブルーチンｂｔｕｎｉｔを呼び出し、そうでない場合には、ｉｓｔａｒｔをｉｓｔａｒｔ２に格納し、ｎｗｉｄｔｈの２分の１をｎｗｉｄｔｈ２に格納し、サブルーチンｂｌｋｔｒｉｄによって３重対角化し、その後、バリア同期を取る。
【００４８】
更に、ｉｓｔａｒｔとｎｗｉｄｔｈ／２を加えたものをｉｓｔａｒｔ３に格納し、ｎｗｉｄｔｈ−ｎｗｉｄｔｈ／２をｎｗｉｄｔｈ３に格納し、ブロックの始まりと終わりをそれぞれ示すｉｓ２、ｉｓ３、ｉｅ２、ｉｅ３、ｉｓ、ｉｅに値を設定し、ｌｅｎ、ｉｐｔｒを設定し、同図の式に従って計算した後、ｕ（ｉｓ：ｉｅ、ｉｓ３：ｉｅ３）に格納して、バリア同期を取り、３重対角化サブルーチンｂｌｋｔｒｉｄを呼び出して処理をしてサブルーチン処理を終わる。
【００４９】
図１４は、３重対角化サブルーチンの内部ルーチンの疑似コードである。
３重対角化の内部サブルーチンｂｔｕｎｉｔにおいては、必要な情報を格納した後、ブロックの始まりｉｐｔｒ２と、幅ｌｅｎ、始まり位置ｉｓ、終わり位置ｉｅを決定し、バリア同期を取る。そして、ｕ（ｉｓ：ｉｅ、ｉ）ｔ＊う（ｉｓ：ｉｅ、ｉ）をｔｍｐに格納し、バリア同期を取る。その後、各値を計算し、それぞれの配列に格納する。ここでｓｕｍは和、ｓｑｒｔは平方根をとることを意味する。そして、最後にバリア同期を取る。
【００５０】
次に、ｖ（ｉｓ：ｉｅ、ｉ）を求めて、バリア同期を取り、ｌｅｎ２、ｉｓｘ、ｉｅｘ、ｕ、ｖを更新して、バリア同期を取る。更に、ｖ（ｉｓ：ｉｅ、ｉ）を更新し、バリア同期を取る。更に、ｖ（ｉｓ：ｉｅ、ｉ）ｔ＊ｕ（ｉｓ：ｉｅ、ｉ）を計算し、ｔｍｐに格納し、バリア同期を取る。
【００５１】
そして、ｂｅｔａに値を設定し、バリア同期を取り、ｖをｂｅｔａを使った演算で更新し、バリア同期を取る。
そして、ｉ＜ｉｂｌｋかつｐｔｒ２＜ｎ−２である時には、ｕ（ｉｓ：ｉｅ、ｉ＋１）を更新し、他の場合には、別の式でｕ（ｉｓ：ｉｅ、ｉ；１：ｉ＋２）を更新し、処理を終了する。なお、このサブルーチンの実行後、スレッドの割当てを解除する。
【００５２】
図１５は、ｕ、ｖから行列の下半分を更新するルーチン及び、対角行列部分の更新をするルーチン及びコピーのルーチンの疑似コードである。
ここで、ｎｂａｓｅはブロックの位置を示すオフセットであり、ｎｗｉｄｔｈは、ブロック幅である。
【００５３】
このサブルーチンｕｐｄａｔｅでは、配列ａ、ｕ、ｖを確保した後、バリア同期を取る。そして、ｂｌｋ、ｎｂａｓｅ２、ｌｅｎ、ｉｓ１、ｉｅ１、ｎｂａｓｅ３、ｉｓｒ、ｉｅｒを設定した後、ａ（ｉｅ１：ｎ、ｉｓ１：ｉｅ１）、及び、ａ（ｉｅｒ＋１：ｎ、ｉｓｒ：ｉｅｒ）をそれぞれ更新し、サブルーチンｔｒｕｐｄａｔｅを２回呼び出し、バリア同期を取って処理を返す。
【００５４】
サブルーチンｔｒｕｐｄａｔｅでは、対角ブロック更新用ブロック幅をｂｌｋ２に格納し、ａ、ｕ、ｖの配列を確保する。
ｉｓ２、ｉｅ２を設定し、Ａ（ｉｓ２：ｉｅ、ｉｓ２：ｉｅ２）を更新し、処理を返す。
【００５５】
サブルーチンｃｏｐｙでは、ｌｅｎ、ｉｓ１、ｌｅｎ１、ｎｂａｓｅ、ｉｓｒ、ｌｅｎｒを設定して、ｂａｎｄｃｐを２回実行し、処理を返す。
図１６は、更新した下三角部分を上三角部分にコピーするルーチンの疑似コードである。
【００５６】
サブルーチンｂａｎｄｃｐでは、ｎｂ、ｗ、ｎｎ、ｌｏｏｐｘを設定し、ｄｏループの中で、ＴＲＬ（ｗ（１：ｎｎｘ、１：ｎｎｘ））にＴＲＬ（ａ（ｉｓ２：ｉｓ２＋ｎｎｘ−１、ｉｓ２：ｉｓ２＋ｎｎｘ））を格納し、ＴＲＵ（ａ（ｉｓ２：ｉｓ２＋ｎｎｘ−１、ｉｓ：ｉｓ＋ｎｎｘ））にＴＲＬ（ｗ（ｉ：ｎｎｘ、１：ｎｎｘ））ｔを格納する。ここで、ＴＲＬは下三角部分、ＴＲＵは上三角部分を示す。
そして、ｗ（１：ｎｎｘ、１：ｎｎｘ）とａ（ｉｓ２：ｉｓ２＋ｎｎｘ、ｉｓ３：ｉｓ３＋ｎｎｘ−１）を更新し、次に、ｗ（１：ｎｙ、１：ｎｘ）とａ（ｉｓ２：ｉｓ２＋ｎｎｘ、ｉｓ３：ｎ）を更新する。
【００５７】
そして、ｄｏループを抜けると、処理を返す。
図１７は、３重対角行列の固有ベクトルを元の行列の固有ベクトルに変換するルーチンの疑似コードである。
【００５８】
ここで、３重対角行列の固有ベクトルは、ｅｖ（１：ｎ、１：ｎｅｖ）に格納されている。また、ａは、３重対角化の出力であり、変換に必要な情報を下三角部分に格納している。
【００５９】
サブルーチンｃｏｎｖｅｖは、配列としてａとｅｖを取る。
スレッドを構成し、処理を並列処理で行う。
スレッド間でバリア同期をとり、ｌｅｎ、ｉｓ、ｉｅ、ｎｅｖｔｈｒｄを設定し、ｃｏｎｖｅｖｔｈｒｄルーチンを呼び出し、リターン後バリア同期を取って処理を終了する。
【００６０】
図１８は、固有ベクトルを変換するルーチンの疑似コードである。
サブルーチンｃｏｎｖｅｖｔｈｒｄでは、ブロック幅をｂｌｋに格納し、配列としてａ、ｅｖ、ｗ、ｗ２を取る。
【００６１】
まず、ｗｉｄｔｈが０より小さければ、処理をせず元のルーチンに戻る。ｎｕｍｂｌｋとｎｆｂｓを設定し、ａｌｐｈａに３重対角化時に対角要素に格納された値の符号を逆にしたもの（−ａ（ｉ、ｉ））を入れる。また、ｘ（１：ｉｗｉｄｔｈ）にｅｖ（ｉ＋１：ｎ、１：ｉｗｉｄｔｈ）ｔ＊ａ（ｉ＋１：ｎ、ｉ）を入力し、これと、ａｌｐｈａとａを用いてｅｖを更新する。更に、次のｄｏ文で、ｉｓ、ｉｅを設定し、ｗ（１：ｂｌｋ、ｉｗｉｄｔｈ）をａ（ｉｓ＋１：ｎ、ｉｓ：ｉｅ）ｔ＊ｅｖ（ｉｓ＋１：ｎ、１：ｉｗｉｄｔｈ）で上書きし、同じくｗ（１：ｂｌｋ、１：ｉｗｉｄｔｈ）をＴＲＬ（ａ（ｉｅ＋１：ｉｓ、ｉｓ：ｉｅ））ｔ＊ｅｖ（ｉｅ＋１：ｉｓ、１：ｉｗｉｄｔｈ）で更新する。ここで、ＴＲＬは、下三角行列部分である。
【００６２】
また、ｗ２の対角要素ベクトルＤＩＡＧ（ｗ２）に、ａ（ｉｓ：ｉｅ、ｉｓ：ｉｅ）の対角要素ベクトルを格納する。
次のｄｏ文でｗ２（ｉ１、ｉ２）をｗ２（ｉ１、ｉ２）＊（ａ（ｉｓ＋ｉ２：ｎ、ｉｓ＋ｉ２−１）ｔ＊ａ（ｉｓ＋ｉ２：ｎ、ｉｓ＋ｉ１−１））で更新する。更に次のｄｏ文で、ｗ２（ｉ１、ｉ２）をｗ２（ｉ１、ｉ２）＋ｗ２（ｉ１、ｉ１＋１：ｉ２−１）＊ｗ２（ｉ１＋１：ｉ２−１、ｉ２）で更新する。
【００６３】
更に、次のｄｏ文で、ｗ２（ｉ１、ｉ２）をｗ２（ｉ１、ｉ２）＊ｗ２（ｉ２、ｉ２）で更新し、ｗ（１：ｂｌｋ、１：ｉｗｉｄｔｈ）、ｅｖ（ｉｓ＋ｎ：ｎ、１：ｉｗｉｄｔｈ）、ｅｖ（ｉｅ＋１：ｉｓ、１：ｉｗｉｄｔｈ）を更新して、処理を元のルーチンに返す。
【００６４】
図１９〜図２９は、疑似コードの処理をフローチャートで示した図である。
図１９は、実対称行列を３重対角化するサブルーチンｔｒｉｄのフローである。ステップＳ１０においては、サブルーチンとしてｓｈａｒｅｄ配列Ａ（ｋ、ｎ）、ｄｉａｇ（ｎ）およびｓｄｉａｇ（ｎ）を入力。ｄｉａｇ、ｓｄｉａｇは計算した３重対角行列の対角要素、副対角要素を出力として返却する。作業域Ｕ（ｎ＋１、ｉｂｌｋ）、ｖ（ｎ＋１、ｉｂｌｋ）をルーチン内部で確保しｓｈａｒｅｄ属性で利用する。ステップＳ１１において、スレッド生成し、各スレッドでローカル域ｎｕｍｔｈｒに総スレッド数ｎｏｔｈｒｄに各スレッドに割り振られたスレッド番号を設定し、各スレッドで以下を設定する。ｉｂｌｋにブロック幅を設定し、ｎｂ＝（ｎ−２＋ｉｂｌｋ−１）／ｉｂｌｋ、ｎｂａｓｅ＝０、ｉ＝１を設定する。ステップＳ１２においては、ｉ＞ｎｂ−１か否かを判断する。ステップＳ１２の判断がＹＥＳの場合には、ステップＳ１９に進む。ステップＳ１２における判断がＮＯの場合には、ステップＳ１３において、ｎｂａｓｅ＝（ｉ−１）×ｉｂｌｋ、ｉｓｔａｒｔ＝１、ｎｗｉｄｔｈ＝ｉｂｌｋを設定する。ステップＳ１４において、サブルーチンｃｏｐｙを呼んで上三角部分に下三角部分をコピーする。ステップＳ１５において、作業領域Ｕにブロック三重対角化の対象領域をコピーする。すなわち、Ｕ（ｎｂａｓｅ＋１：ｎ、１：ｉｂｌｋ）←Ａ（ｎｂａｓｅ＋１：ｎ、ｎｂａｓｅ＋１：ｎｂａｓｅ＋ｉｂｌｋ）を実行する。ステップＳ１６において、サブルーチンｂｌｋｔｒｉｄを呼んでＵにコピーした部分の三重対角化を行う。（ｉｓｔａｒｔ＝１、ブロック幅はｉｂｌｋを受け渡す）。ステップＳ１７において、三重対角化されたものを配列Ａに戻す。すなわち、Ａ（ｎｂａｓｅ＋１：ｎ，ｎｂａｓｅ＋１：ｎｂａｓｅ＋ｉｂｌｋ）←Ｕ（ｎｂａｓｅ＋１：ｎ、１：ｉｂｌｋ）を行う。ステップＳ１８において、サブルーチンｕｐｄａｔｅを呼んで、Ａ（ｎｂａｓｅ＋ｉｂｌｋ：ｎ、ｎｂａｓｅ＋ｉｂｌｋ：ｎ）の下三角部分を更新し、ステップＳ１２に戻る。
【００６５】
ステップＳ１９においては、ｎｂａｓｅ＝（ｎｂ−１）×ｉｂｌｋ、ｉｓｔａｒｔ＝１、ｉｂｌｋ２＝ｎ−ｎｂａｓｅと設定する。ステップＳ２０においては、作業領域Ｕにブロック三重対角化の対象領域をコピーする。すなわち、Ｕ（ｎｂａｓｅ＋１：ｎ，１：ｎｗｉｄｔｈ）←Ａ（ｎｂａｓｅ＋１：ｎ、ｎｂａｓｅ＋１：ｎ）を行う。ステップＳ２１において、サブルーチンｂｌｋｔｒｉｄを呼んでＵにコピーした部分の三重対角化を行う（ｉｓｔａｒｔ＝１、ブロック幅はｉｂｌｋ２を受け渡す）。ステップＳ２２において、三重対角化されたものを配列Ａに戻す。すなわち、Ａ（ｎｂａｓｅ＋１：ｎ、ｎｂａｓｅ＋１：ｎ）←Ｕ（ｎｂａｓｅ＋１：ｎ、１：ｎｗｉｄｔｈ）を行う。ステップＳ２３においては、並列処理のために生成したスレッドを消して、サブルーチンを終了する。
【００６６】
図２０は、サブルーチンｂｌｋｔｒｉｄのフローである。このサブルーチンは再帰的プログラムとなっている。
このサブルーチンは以下のステートメントによって呼び出される。
ｓｕｂｒｏｕｔｉｎｅ　ｂｌｋｔｒｉｄ（Ａ，ｋ，ｎ，ｄｉｇ，ｓｄｉｇ，ｎｂａｓｅ，ｉｓｔａｒｔ，ｎｗｉｄｔｈ、Ｕ、Ｖ、ｎｏｔｈｒｄ、ｎｕｍｔｈｒｄ）
ここで、ｎｂａｓｅはブロックの位置を示すオフセット。ｉｓｔａｒｔは再帰呼び出しで対称となる縮小されたブロックのブロック内でのオフセットで最初は１で再帰的に呼び出されるとき対象ブロックの位置を示す。ｎｗｉｄｔｈはブロック幅を表す。ステップＳ２５においては、ｎｗｉｄｔｈ＜１０か否かを判断する。ステップＳ２５の判断がＮＯの場合には、ステップＳ２７に進む。ステップＳ２５の判断がＹＥＳの場合には、ステップＳ２６において、サブルーチンｂｔｕｎｉｔを呼んで３重対角化を行い、サブルーチンを抜ける。ステップＳ２７では、再帰呼び出しのために対象となる更新位置とブロック幅を変えて、呼び出すｉｓｔａｒｔ２＝ｉｓｔａｒｔ、ｎｗｉｄｔｈ２＝ｎｗｉｄｔｈ／２を設定し受け渡す。縮小されたブロックの開始位置、ブロック幅を受け渡す。ステップＳ２８においては、サブルーチンｂｌｋｔｒｉｄを再帰的に呼び出す。ステップＳ２９においては、スレッド間でバリア同期を取る。ステップＳ３０においては、更新で各スレッドが分担する始点（ｉｓ２、ｉｓ３）、終点（ｉｅ２、ｉｅ３）を計算する。すなわち、ｉｓｔａｒｔ３＝ｉｓｔａｒｔ＋ｎｗｉｄｔｈ／２、ｎｗｉｄｔｈ３＝ｎｗｉｄｔｈ−ｎｗｉｄｔｈ／２、ｉｓ２＝ｉｓｔａｒｔ２、ｉｅ２＝ｉｓｔａｒｔ＋ｎｗｉｄｔｈ２−１、ｉｓ３＝ｉｓｔａｒｔ３、ｉｅ３＝ｉｓｔａｒｔ３＋ｎｗｉｄｔｈ３−１、ｉｐｔｒ＝ｎｂａｓｅ＋ｉｓｔａｒｔ３、ｌｅｎ＝（ｎ−ｉｐｔｒ＋ｎｕｍｔｈｒｄ−１）／ｎｕｍｔｈｒｄ、ｉｓ＝ｉｐｔｒ＋（ｎｏｔｈｒｄ−１）×ｌｅｎ＋１、ｉｅ＝ｍｉｎ（ｎ，ｉｐｔｒ＋ｎｏｔｈｒｄ×ｌｅｎ）を計算する。ステップＳ３１においては、Ｕ（ｉｓ：ｉｅ，ｉｓ３：ｉｅ３）＝Ｕ（ｉｓ：ｉｅ，ｉｓ３：ｉｅ３）−Ｕ（ｉｓ：ｉｅ，ｉｓ２：ｉｅ２）×Ｗ（ｉｓ３：ｉｅ３，ｉｓ２：ｉｅ２）ｔ−Ｗ（ｉｓ：ｉｅ，ｉｓ２：ｉｅ２）×Ｕ（ｉｓ３：ｉｅ３，ｉｓ２：ｉｅ２）ｔを計算し、ステップＳ３２において、スレッド間でバリア同期を取り、ステップＳ３３においてサブルーチンｂｌｋｔｒｉｄを再帰的に呼び出して、サブルーチンを抜ける。
【００６７】
図２１及び図２２は、サブルーチンｂｌｋｔｒｉｄの内部ルーチンであるサブルーチンｂｔｕｎｉｔのフローである。
ステップＳ３５において、ｔｍｐ（ｎｕｍｔｈｒｄ）、ｓｉｇｍａ、ａｌｐｈａをｓｈａｒｅｄ属性で割り付ける。ステップＳ３６において、ｎｂａｓｅ＋ｉｓｔａｒｔ＞ｎ−２か否かを判断する。ステップＳ３６の判断がＹＥＳの場合には、サブルーチンを抜ける。ステップＳ３６の判断がＮＯの場合にはステップＳ３８に進む。ステップＳ３８では、ｉ＝ｉｓｔａｒｔとし、ステップＳ３９において、ｉ＜＝ｉｓｔａｒｔ−１＋ｎｗｉｄｔｈであるか否かを判断する。ステップＳ３９の判断がＮＯの場合には、サブルーチンを抜ける。ステップＳ３９の判断がＹＥＳの場合には、ステップＳ４０において、各スレッドで分担する始点ｉｓ、終点ｉｅを計算する。ｉｐｔｒ２＝ｎｂａｓｅ＋ｉ、ｌｅｎ＝（ｎ−ｉｐｔｒ２＋ｎｕｍｔｈｒｄ−１）／ｎｕｍｔｈｒｄ、ｉｓ＝ｉｐｔｒ２＋（ｎｏｔｈｒｄ−１）×ｌｅｎ＋１、ｉｅ＝ｍｉｎ（ｎ，ｉｐｔｒ２＋ｎｏｔｈｒｄ×ｌｅｎ）を計算する。ステップＳ４１では、バリア同期を取り、ステップＳ４２では、ｔｍｐ（ｎｏｔｈｒｄ）＝Ｕ（ｉｓ：ｉｅ，ｉ）ｔ×Ｕ（ｉｓ：ｉｅ，ｉ）を計算し、ステップＳ４３で、バリア同期を取る。ステップＳ４４において、ｎｏｔｈｒｄ＝１か否かを判断する。ステップＳ４４の判断がＮＯの場合には、ステップＳ４６に進む。ステップＳ４４の判断がＹＥＳの場合には、ステップＳ４５において、各スレッドで部分計算した和の平方根をとり、三重対角化のための計算を行う（ハウスホルダーベクトルの作成）。ｓｉｇｍａ＝ｓｑｒｔ（ｓｕｍ（ｔｍｐ（１：ｎｕｍｔｈｒｄ）））、ここで、ＳＵＭは和、ｓｑｒｔは平方根である。ｄｉａｇ（ｉｐｔｒ２）＝ｕ（ｉｐｔｒ２，ｉ）、ｓｄｉａｇ（ｉｐｔｒ２）＝−ｓｉｇｍａ、Ｕ（ｎｂａｓｅ＋ｉ＋１、ｉ）＝Ｕ（ｎｂａｓｅ＋ｉ＋１，ｉ）＋ｓｉｇｎ（ｕ（ｎｂａｓｅ＋ｉ＋１，ｉ）×ｓｉｇｍａ、ａｌｐｈａ＝１．０／（ｓｉｇｍａ×ｕ（ｎｂａｓｅ＋ｉ＋１，ｉ）、Ｕ（ｉｐｔｒ２，ｉ）＝ａｌｐｈａを計算して、ステップＳ４６に進む。ステップＳ４６では、バリア同期を取り、ステップＳ４７では、ｉｐｔｒ３＝ｉｐｔｒ２＋１、ステップＳ４８では、Ｖ（ｉｓ：ｉｅ，ｉ）＝Ａ（ｉｐｔｒ３：ｎ、ｉｐｔｒ２＋ｉｓ：ｉｐｔｒ２＋ｉｅ）ｔＵ（ｐｔｒ３：ｎ，ｉ）を演算し、ステップＳ４９では、バリア同期を取る。
【００６８】
ステップＳ５０では、Ｖ（ｉｓ：ｉｅ，ｉ）＝ａｌｐｈａ×Ｖ（ｉｓ：ｉｅ，ｉ）−Ｖ（ｉｓ：ｉｅ，１：ｉ−１）×（Ｕ（ｉｐｔｒ３：ｎ，１：ｉ−１）ｔ×Ｕ（ｉｐｔｒ３：ｎ、ｉ））−Ｕ（ｉｓ：ｉｅ，１：ｉ−１）×（Ｖ（ｉｐｔｒ３：ｎ，１：ｉ−１）ｔ×Ｕ（ｉｐｔｒ３：ｎ，ｉ））を計算し、ステップＳ５１において、バリア同期を取る。ステップＳ５２においては、ｔｍｐ（ｎｏｔｈｒｄ）＝Ｖ（ｉｓ：ｉｅ，ｉ）ｔ×Ｕ（ｉｓ：ｉｅ，ｉ）を計算し、ステップＳ５３において、バリア同期を取る。ステップＳ５４においては、ｎｏｔｈｒｄ＝１か否かを判断する。ステップＳ５４の判断がＮＯの場合には、ステップＳ５６に進み、ステップＳ５４の判断がＹＥＳの場合には、ステップＳ５５に進む。ステップＳ５５においては、ｂｅｔａ＝０．５×ａｌｐｈａ×ｓｕｍ（ｔｍｐ（１：ｎｕｍｔｈｒｄ））を計算する。ただし、ｓｕｍはベクトルの和をとる記号である。ステップＳ５６では、バリア同期を取る。ステップＳ５７では、Ｖ（ｉｓ：ｉｅ、ｉ）＝Ｖ（ｉｓ：ｉｅ，ｉ）−ｂｅｔａ×Ｕ（ｉｓ：ｉｅ，ｉ）を計算し、ステップＳ５８において、バリア同期を取る。ステップＳ５９においては、ｐｔｒ２＜ｎ−２か否かを判断し、判断がＹＥＳの場合には、ステップＳ６０において、Ｕ（ｉｓ：ｉｅ，ｉ＋１）＝Ｕ（ｉｓ：ｉｅ，ｉ＋１）−Ｕ（ｉｓ：ｉｅ，ｉｓｔａｒｔ：ｉ）×Ｖ（ｉ＋１，ｉｓｔａｒｔ：１）ｔ−Ｖ（ｉｓ：ｉｅ，ｉｓｔａｒｔ：ｉ）×Ｕ（ｎ＋１，ｉｓｔａｒｔ：１）ｔを計算し、ステップＳ３９に戻る。ステップＳ５９の判断がＮＯの場合には、ステップＳ６１において、Ｕ（ｉｓ：ｉｅ，ｉ＋１：ｉ＋２）＝Ｕ（ｉｓ：ｉｅ，ｉ＋１：Ｉ＋２）−Ｕ（ｉｓ：ｉｅ，ｉｓｔａｒｔ：ｉ）×Ｖ（ｉ＋１：ｎ，ｉｓｔａｒｔ：ｉ）ｔ−Ｖ（ｉｓ：ｉｅ，ｉｓｔａｒｔ：ｉ）×Ｕ（ｎ＋１：ｎ，ｉｓｔａｒｔ：ｉ）ｔを計算して、サブルーチンを抜ける。
【００６９】
図２３は、サブルーチンｕｐｄａｔｅのフローである。
ステップＳ６５において、バリア同期を取り、ステップＳ６６において、各スレッドでペアを作り、更新を分担する始点、終点を決める。すなわち、ｎｂａｓｅ２＝ｎｂａｓｅ＋ｉｂｌｋ、ｌｅｎ＝（ｎ−ｎｂａｓｅ２＋２×ｎｕｍｔｈｒｄ−１）／（２×ｎｕｍｔｈｒｄ）、ｉｓ１＝ｎｂａｓｅ２＋（ｎｏｔｈｒｄ−１）ｌｅｎ＋１、ｉｅ１＝ｍｉｎ（ｎ，ｎｂａｓｅ２＋ｎｏｔｈｒｄ×ｌｅｎ）、ｎｂａｓｅ３＝ｎｂａｓｅ２＋２×ｎｕｍｔｈｒｄ×ｌｅｎ、ｉｓｒ＝ｎｂａｓｅ３−ｎｏｔｈｒｄ×ｌｅｎ＋１、ｉｅｒ＝ｍｉｎ（ｎ，ｉｓｒ＋ｌｅｎ−１）を計算する。ステップＳ６７では、Ａ（ｉｅ１＋１：ｎ，ｉｓ１：ｉｅ１）＝Ａ（ｉｅ１＋１：ｎ，ｉｓ１＋１：ｎ，ｉｓ１：ｉｅ１）−Ｗ（ｉｅ１＋１：ｎ，１：ｂｌｋ）×Ｕ（ｉｓ１：ｉｅ１，１：ｂｌｋ）ｔ−Ｕ（ｉｅ１＋１：ｎ，１：ｂｌｋ）×Ｗ（ｉｓ１：ｉｅ１，１：ｂｌｋ）ｔ、Ａ（ｉｅｒ＋１：ｎ，ｉｓｒ：ｉｅｒ）＝Ａ（ｉｅｒ＋１：ｎ，ｉｓｒ：ｉｅｒ）−Ｗ（ｉｅｒ＋１：ｎ，１：ｂｌｋ）×Ｕ（ｉｓｒ：ｉｅｒ，１：ｂｌｋ）ｔ−Ｕ（ｉｅｒ＋１：ｎ，１：ｂｌｋ）×Ｗ（ｉｓｒ：ｉｅｒ，１：ｂｌｋ）ｔを計算し、ステップＳ６８で、サブルーチンｔｒｕｐｄａｔｅを呼んで左半分の対角行列部分を更新する。ｉｓ１、ｉｅ１、Ａ、Ｗ、Ｕを渡す。ステップＳ６９において、サブルーチンｔｒｕｐｄａｔｅを呼んで右半分の対角行列部分を更新する。ｉｓｒ、ｉｅｒ、Ａ、Ｗ、Ｕを渡す。ステップＳ７０において、バリア同期を取って、サブルーチンを抜ける。
【００７０】
図２４は、サブルーチンｔｒｕｐｄａｔｅ（対角行列部分の更新）のフローである。
更新の始点ｉｓ、と終点ｉｅを入力とする。これは、サブルーチンが呼び出される前に対角ブロックの下の長方形の部分の更新に使われた始点と終点である。
【００７１】
ステップＳ７５においてｂｌｋ２に対角ブロック更新用ブロック幅を設定し、ｉ＝ｉｓを設定する。ステップＳ７６において、ｉ＞ｉｅ−１か否かを判断する。ステップＳ７６における判断がＹＥＳの場合には、サブルーチンを抜ける。ステップＳ７６における判断がＮＯの場合には、ステップＳ７７において、各スレッドでの更新の始点、終点を決める。すなわち、ｉｓ２＝ｉ、ｉｅ２＝ｍｉｎ（ｉ＋ｂｌｋ２−１，ｉｅ−１）、Ａ（ｉｓ２：ｉｅ−１，ｉｓ２，ｉｅ２）＝Ａ（ｉｓ２：ｉｅ−１，ｉｓ２，ｉｅ２）−Ｕ（ｉｓ２：ｉｅ−１，１：ｂｌｋ）×Ｗ（ｉｓ２：ｉｅ２，１：ｂｌｋ）ｔ−Ｗ（ｉｓ２：ｉｅ−１，１：ｂｌｋ）×Ｕ（ｉｓ２：ｉｅ２，１：ｂｌｋ）ｔを計算し、ステップＳ７８において、ｉ＝ｉ＋ｂｌｋ２と設定し、ステップＳ７６に戻る。
【００７２】
図２５は、サブルーチンｃｏｐｙのフローである。
ステップＳ８０において、コピーを各スレッドで対を作って並列実行するときの始点、と幅を計算する。すなわち、ｌｅｎ＝（ｎ−ｎｂａｓｅ＋２×ｎｕｍｔｈｒｄ−１）／（２×ｎｕｍｔｈｒｄ）、ｉｓ１＝ｎｂａｓｅ＋（ｎｏｔｈｒｄ−１）×ｌｅｎ＋１、ｌｅｎ１＝ｍａｘ（０，ｍｉｎ（ｎ−ｉｓ１＋１，ｌｅｎ））、ｎｂａｓｅ３＝ｎｂａｓｅ＋２×ｎｕｍｔｈｒｄ×ｌｅｎ、ｉｓｒ＝ｎｂａｓｅ３−ｎｏｔｈｒｄ×ｌｅｎ＋１、ｌｅｎｒ＝ｍａｘ（０，ｍｉｎ（ｎ−ｉｓｒ＋１，ｌｅｎ））を計算する。ステップＳ８１において、サブルーチンｂａｎｄｃｐを呼び出す。ペアの左側の始点ｉｓ１、幅ｌｅｎ１できまる領域をコピーする。ステップＳ８２において、サブルーチンｂａｎｄｃｐを呼び出す。ペアの右側の始点ｉｓｒ、幅ｌｅｎｒできまる領域をコピーする。
【００７３】
図２６は、サブルーチンｂａｎｄｃｐのフローである。
このルーチンは、小さな作業域ＷＸを利用してキャッシュ上での転置を行いながらコピーを行う。作業域をＷＸ（ｎｂ、ｎｂ）とし、始点をｉｓ、幅をｌｅｎで受ける。
【００７４】
ステップＳ８５において、ｎｎ＝ｍｉｎ（ｎｂ，ｌｅｎ）、ｌｏｏｐｘ＝（ｌｅｎ＋ｎｎ−１）／ｎｎ、ｊ＝１を演算する。ステップＳ８６において、ｊ＞ｌｏｏｐｘか否かを判断する。ステップＳ８６の判断がＹＥＳの時は、サブルーチンを終了する。ステップＳ８６の判断がＮＯの場合には、ステップＳ８７において、ＷＸにコピーする対角ブロックの大きさｎｎｘとそのオフセットｉｐを決める。ｉｐ＝ｉｓ＋（ｊ−１）×ｎｎ、ｎ１＝ｌｅｎ−（ｊ−１）×ｎｎ、ｎｎｘ＝ｍｉｎ（ｎｎ，ｎ１）、ｌｅｎ２＝ｎ−ｉｐ−ｎｎｘ＋１、ｌｏｏｐｙ＝（ｌｅｎ２＋ｎｎ−１）／ｎｎ、ＴＲＬ（ＷＸ（１：ｎｎｘ，１：ｎｎｘ））＝ＴＲＬ（Ａ（ｉｐ：ｉｐ＋ｎｎｘ−１，ｉｐ：ｉｐ＋ｎｎｘ−１））、ＴＲＵ（Ａ（ｉｐ：ｉｐ＋ｎｎｘ−１，ｉｐ：ｉｐ＋ｎｎｘ−１））＝ＴＲＬ（ＷＸ（１：ｎｎｘ，ｉ：ｎｎｘ））、ｉ＝１，ｉｓ２＝ｉｐ，ｉｓ３＝ｉｐ＋ｎｎｘを計算する。ただし、ＴＲＵは上三角部分、ＴＲＬは下三角部分。
【００７５】
ステップＳ８８では、ｉ＞ｌｏｏｐｙ−１か否かを判断する。ステップＳ８８の判断がＮＯの場合には、ステップＳ８９において、ｎｎ×ｎｎｘの部分を転置してコピーする。すなわち、ＷＸ（１：ｎｎ、１：ｎｎｘ）＝Ａ（ｉｓ３：ｉｓ３＋ｎｎ−１、ｉｓ２：ｉｓ２＋ｎｎｘ−１）、Ａ（ｉｓ２：ｉｓ２＋ｎｎｘ−１，ｉｓ３：ｉｓ３＋ｎｎ−１）＝ＷＸ（１，ｎｎ：１，ｎｎｘ）ｔ、ｉｓ３＝ｉｓ３＋ｎｎを計算して、ステップＳ８８に戻る。ステップＳ８８の判断がＹＥＳの場合には、ステップＳ９０で、最後の部分をコピーする。すなわち、ｎｎ＝ｎ−ｉｓ３＋１、ＷＸ（１：ｎｎ，１：ｎｘ）＝Ａ（ｉｓ３：ｎ，ｉｓ２：ｉｓ２＋ｎｎｘ−１）、Ａ（ｉｓ２：ｉｓ２＋ｎｎｘ−１，ｉｓ３：ｎ）＝ＷＸ（１：ｎｎ，１：ｎｘ）を演算して、ステップＳ８６に戻る。
【００７６】
図２７は、サブルーチンｃｏｎｖｅｖのフローである。
このルーチンでは、求める固有ベクトルの数ｎｅｖ、ハウスホルダーベクトルはａの下半分に格納されている。３重対角行列の固有ベクトルがｅｖ（ｋ，ｎｅｖ）に格納されている。
【００７７】
ステップＳ９５において、スレッドを生成する。各スレッドのローカル域のｎｕｍｔｈｒ及びｎｏｔｈｒｄに総スレッド数とスレッド番号（１〜ｎｕｍｔｈｒｄ）を設定する。ステップＳ９６において、バリア同期を取る。ステップＳ９７において、各スレッドで分担計算する始点、終点を決める。すなわち、ｌｅｎ＝（ｎｅｖ＋ｎｕｍｔｈｒｄ−１）／ｎｕｍｔｈｒｄ、ｉｓ＝（ｎｏｔｈｒｄ−１）×ｌｅｎ＋１、ｉｅ＝ｍｉｎ（ｎｅｖ，ｎｏｔｈｒｄ×ｌｅｎ）、ｗｉｄｔｈ＝ｉｅ−ｉｓ＋１を計算する。ステップＳ９８において、サブルーチンｃｏｎｖｅｖｔｈｒｄを呼び出し、三重対角行列からもとの行列の固有ベクトルに変換する。各スレッドが担当する固有ベクトルの格納されている領域と固有ベクトルの数ｗｉｄｔｈを渡す。ステップＳ９９において、バリア同期をとり、ステップＳ１００において、生成したスレッドを消滅させ、サブルーチンを抜ける。
【００７８】
図２８及び図２９は、サブルーチンｃｏｎｖｅｖｔｈｒｄのフローである。
このルーチンは、各スレッドごとに分担した三重対角行列の固有ベクトルを変換して、もとの行列の固有ベクトルにする。ハウスホルダ変換を復元するベクトルと係数は配列Ａに格納されている。
【００７９】
ステップＳ１１０において、ｂｌｋにブロック幅を設定する。ブロック幅は８０程度である。ステップＳ１１１では、ｉｗｉｄｔｈ＜０か否かを判断する。ステップＳ１１１の判断がＹＥＳの場合には、サブルーチンを抜ける。ステップＳ１１１の判断がＮＯの場合には、ステップＳ１１２に進む。ステップＳ１１２では、以下のループで変換の最初のブロックは、（１＋αｕｕｔ）を順次計算して求める。まず、ｎｕｍｂｌｋ＝（ｎ−２＋ｂｌｋ−１）／ｂｌｋ、ｎｆｂｓ＝ｎ−２−ｂｌｋ×（ｎｕｍｂｌｋ−１）を求める。ステップＳ１１３において、ｉ＜ｎ−２−ｎｆｂｓ＋１か否かを判断する。ステップＳ１１３の判断がＹＥＳの場合には、ステップＳ１１４に進み、ａｌｐｈａ＝−ａ（ｉ，ｉ）、ｘ（１：ｉｗｉｄｔｈ）＝ａ（ｉ＋１：ｎ，ｉ）ｔ×ｅｖ（ｉ＋１：ｎ，１：ｗｉｄｔｈ）、ｅｖ（ｉ＋１：ｎ，１：ｗｉｄｔｈ）＝ｅｖ（ｉ＋１：ｎ，１：ｗｉｄｔｈ）＋ａｌｐｈａ×ａ（ｉ＋１：ｎ，ｉ）×ｘ（１：ｉｗｉｄｔｈ）ｔを演算し、ステップＳ１１３に戻る。ステップＳ１１３の判断がＮＯの場合には、ステップＳ１１５において、ｉ＝１に設定し、ステップＳ１１６において、ｉ＞ｎｕｍｂｌｋ−１か否かを判断する。ステップＳ１１６の判断がＮＯの場合には、サブルーチンを抜ける。ステップＳ１１６の判断がＹＥＳの場合には、ステップＳ１１７において、ブロック形式の（１＋ＵＢＵｔ）のＵｔ×ＥＶをＵｔの左端の上三角行列部分と右側の長方形の部分に分けて計算する。すなわち、ｉｓ＝ｎ−２−（ｎｆｎｓ＋ｉ×ｂｌｋ）＋１、ｉｅ＝ｉｅ＋ｂｌｋ−１、Ｗ（１：ｂｌｋ，ｉｗｉｄｔｈ）＝ａ（ｉｅ＋１：ｎ，ｉｓ：ｉｅ）ｔ×ｅｖ（ｉｅ＋１：ｎ，１：ｉｗｉｄｔｈ）、Ｗ（１：ｂｌｋ−１，１：ｉｗｉｄｔｈ）＝ｗ（１：ｂｌｋ−１，１：ｉｗｉｄｔｈ）＋ＴＲＬ（ａ（ｉｓ＋１：ｉｅ，ｉｓ：ｉｅ−１））ｔ×ｅｖ（ｉｓ＋１：ｉｅ，１：ｉｗｉｄｔｈ）を計算する。次に、ブロック形式の（１＋ＵＢＵｔ）のＢを求める。ｄｉａｇ（ｗ２）＝−ｄｉａｇ（ａ（ｉｓ：ｉｓ＋ｂｌｋ−１，ｉｓ：ｉｓ＋ｂｌｋ−１））、ｉ２＝ｂｌｋ、を計算する。ｗ２の部分に対応する係数αを格納する。なお、ＴＲＬ（ｗ２）はｗ２の下三角行列、ｄｉａｇ（ｘ）はｘの対角要素部分である。
【００８０】
ステップＳ１１８においては、ｉ２＜１か否かを判断する。ステップＳ１１８の判断がＮＯの場合には、ステップＳ１１９において、ｗ２の上三角部分にハウスホルダベクトルの内積×αを格納し、ｉ１＝ｉ２−１と設定する。ステップＳ１２０では、ｉ１＜１か否かを判断し、ＮＯの場合には、ステップＳ１２１において、ｗ２（ｉ１，ｉ２）＝ｗ２（ｉ１，ｉ１）×（ａ（ｉｓ＋ｉ２：ｎ，ｉｓ＋ｉ２−１）ｔ×ａ（ｉｓ＋ｉ２：ｎ，ｉｓ＋ｉ１−１））、ｉ１＝ｉ１−１を計算し、ステップＳ１２０に戻る。ステップＳ１２０の判断がＹＥＳの場合には、ステップＳ１２２において、ｉ２＝ｉ２−１と設定し、ステップＳ１１８に戻る。ステップＳ１１８の判断がＹＥＳの場合には、ステップＳ１２３において、ｉ１＝ｂｌｋ−２とし、以下、２重ループで展開係数を計算する。三角行列の上辺を右から左に決定し、積み上げるように計算する。これは左からハウスホルダ変換をかけたときの展開を足して係数を決めることに相当する。ステップＳ１２４において、ｉ１＜１か否かを判断する。ステップＳ１２４の判断がＮＯの場合には、ステップＳ１２５においては、ｉ２＝ｂｌｋと設定し、ステップＳ１２６においては、ｉ２＜ｉ１＋１か否かを判断する。ステップＳ１２６の判断がＮＯの場合には、ステップＳ１２７において、上辺の要素を左から右に決めていく。一つ前の係数を利用する。すなわち、ｗ２（ｉ１，ｉ２）＝ｗ２（ｉ１，ｉ２）＋ｗ２（ｉ１，ｉ１＋１：ｉ２−１）×ｗ２（ｉ１＋１：ｉ２−１，ｉ２）、ｉ２＝ｉ２−１を計算し、ステップＳ１２６に戻る。ステップＳ１２６の判断がＹＥＳの場合には、ステップＳ１２８において、ｉ１＝ｉ１−１と設定し、ステップＳ１２４に戻る。ステップＳ１２４の判断がＹＥＳの場合には、ステップＳ１２９に進み、ｉ２＝ｂｌｋを設定し、ステップＳ１３０において、ｉ２＜１か否かを判断する。ステップＳ１３０の判断がＮＯの場合には、ステップＳ１３１において、以下のループで不足していた係数αをかける。まず、ｉ１＝ｉ２−１と設定し、ステップＳ１３２において、ｉ１＜１か否かを判断する。ステップＳ１３２における判断がＮＯの場合には、ステップＳ１３３において、ｗ２（ｉ１，ｉ２）＝ｗ２（ｉ１，ｉ２）×ｗ２（ｉ２，ｉ２）、ｉ１＝ｉ１−１を計算し、ステップＳ１３２に戻る。ステップＳ１３２の判断がＹＥＳの場合には、ステップＳ１３４において、ｉ２＝ｉ２−１と設定し、ステップＳ１３０に戻る。ステップＳ１３０の判断がＹＥＳの場合、ステップＳ１３５において、ＢＵｔを計算してＷに格納する。Ｗ（１：ｂｌｋ，１：ｉｗｉｄｔｈ）＝ＴＲＵ（ｗ２）×Ｗ（１：ｂｌｋ，１：ｉｗｉｄｔｈ）を計算する。そして、Ｕの上部の三角行列と下部の長方形の部分とＷに格納したＢＵｔを使って（１＋ＵＢＵｔ）×ＥＶを計算する。すなわち、ｅｖ（ｉｅ＋１：ｎ，１：ｗｉｄｔｈ）＝ｅｖ（ｉｅ＋１：ｎ，１：ｗｉｄｔｈ）＋ａ（ｉｅ＋１：ｎ，ｉｓ：ｉｅ）×Ｗ（１：ｂｌｋ，１：ｗｉｄｔｈ）、ｅｖ（ｉｓ＋１：ｉｅ，１：ｗｉｄｔｈ）＝ｅｖ（ｉｓ＋１：ｉｅ，１：ｗｉｄｔｈ）＋ＴＲＬ（ａ（ｉｓ＋１：ｉｅ，ｉｓ＋１：ｉｅ））×Ｗ（１：ｂｌｋ−１，１：ｗｉｄｔｈ）を計算し、ステップＳ１１５に戻る。
【００８１】
本発明によれば、共有メモリ型スカラ並列計算機を用いて高性能かつスケーラビリティのある固有値・固有ベクトルの並列演算方法を提供できる。
本発明の実施形態によれば、特に、固有ベクトルの変換に関する計算は従来の方法に比べて１０倍程度スカラマシンで高速。これらのアルゴリズムを使った実対称行列の固有値・固有ベクトルは、３重対角行列の固有値・固有ベクトルを同じスツルムの方法と逆反復法で求める方法を利用したもので、ＳＵＮの数値計算ライブラリＳＵＮ　Ｐｅｒｆｏｒｍａｎｃｅ　ｌｉｂｒａｒｙの機能に比べて、７個のＣＰＵで６．７倍高速。ＳＵＮの他のルーチン（機能が劣る：固有値・固有ベクトルを選択的には求められない）３重対角行列の固有値・固有ベクトルの計算をｄｉｖｉｄｅ　＆　ｃｏｎｑｕｅｒ法で行うものと比べても２．３倍高速。
【００８２】
これらのアルゴリズムを使ったエルミート行列の固有値・固有ベクトルは、３重対角行列の固有値・固有ベクトルを同じスツルムの方法と逆反復法で求める方法を利用したもので、ＳＵＮの数値計算ライブラリＳＵＮ　Ｐｅｒｆｏｒｍａｎｃｅ　ｌｉｂｒａｒｙの機能に比べて７個のＣＰＵで４．８倍高速。ＳＵＮの他のルーチン（機能が劣る：固有値を選択的には求められない）３重対角行列の固有値・固有ベクトルの計算をｄｉｖｉｄｅ　＆　ｃｏｎｑｕｅｒ法で行うものと比べても３．８倍高速。
【００８３】
なお、行列計算の基本的アルゴリズムに関しては、以下の教科書を参照されたい。
Ｇ．　Ｈ．　Ｇｏｌｕｂ　ａｎｄ　Ｃ．　Ｆ．　Ｖａｎ　Ｌｏａｎ　”Ｍａｔｒｉｘ　Ｃｏｍｐｕｔａｔｉｏｎｓ”　ｔｈｉｒｄ　ｅｄｉｔｉｏｎ　Ｔｈｅ　Ｊｏｈｎｓ　Ｈｏｐｋｉｎｓ　Ｕｎｉｖｅｒｓｉｔｙ　Ｐｒｅｓｓ，　１９９６
また、３重対角化の並列計算については、以下の文献を参照されたい。
Ｊ．　Ｃｈｏｉ，　Ｊ．　Ｊ．　Ｄｏｎｇａｒｒａ，　ａｎｄ　Ｄ．　Ｗ．　Ｗａｌｋｅｒ，　”ＴＨＥ　ＤＥＳＩＧＮ　ＯＦ　Ａ　ＰＡＲＡＬＬＥＬ　ＤＥＮＳＥ　ＬＩＮＥＡＲ　ＡＬＧＥＢＲＡ　ＳＯＦＴＷＡＲＥ　ＬＩＢＲＡＲＹ：　ＲＥＤＵＣＴＩＯＮ　ＴＯ　ＨＥＳＳＥＮＢＥＲＧ，　ＴＲＩＤＩＡＧＯＮＡＬ，　ＡＮＤ　ＢＩＤＩＡＧＯＮＡＬ　ＦＯＲＭ”，　Ｅｎｇｉｎｅｅｒｉｎｇ　Ｐｈｙｓｉｃｓ　ａｎｄ　Ｍａｔｈｅｍａｔｉｃｓ　Ｄｉｖｉｓｉｏｎ，　Ｍａｔｈｅｍａｔｉｃａｌ　Ｓｃｉｅｎｃｅｓ　Ｓｅｃｔｉｏｎ，　ｐｒｅｐａｒｅｄ　ｂｙ　ｔｈｅ　Ｏａｋ　Ｒｉｄｇｅ　Ｎａｔｉｏｎａｌ　Ｌａｂｏｒａｔｏｒｙ　ｍａｎａｇｅｄ　ｂｙ　Ｍａｒｔｉｎ　Ｍａｒｉｅｔｔａ　Ｅｎｅｒｇｙ　Ｓｙｓｔｅｍ，　Ｉｎｃ．，　ｆｏｒ　ｔｈｅ　Ｕ．Ｓ．　ＤＥＰＡＲＴＭＥＮＴ　ＯＦ　ＥＮＥＲＧＹ　ｕｎｄｅｒ　Ｃｏｎｔｒａｃｔ　Ｎｏ．　ＤＥ−ＡＣ０５−８４ＯＲ２１４００，　ＯＲＮＬ／ＴＭ−１２４７２．
（付記１）共有メモリ型スカラ並列計算機用固有値問題の並列処理方法を実現させるプログラムであって、
演算すべき実対称行列あるいはエルミート行列をブロック化し、ブロック化した部分をメモリの作業域にコピーして、ブロック間の積を利用して３重対角化するステップと、
３重対角化された行列から固有値と固有ベクトルを求めるステップと、
該３重対角化された行列から求められた固有ベクトルを、所定のブロック幅の行列演算の並列演算となるように変形されたハウスホルダー変換によって変換し、元の行列の固有ベクトルを求めるステップと、
を備えることを特徴とする並列処理方法を共有メモリ型スカラ並列計算機に実現させるプログラム。
【００８４】
（付記２）前記３重対角化ステップにおいて、ブロック化された部分の更新を再帰的プログラムで行うことを特徴とする付記１に記載のプログラム。
（付記３）前記３重対角化ステップにおいて、ブロック化された部分は、キャッシュメモリのページを跨らないように小さなブロックに分割して格納し、各プロセッサで並列に演算を行うことによって処理されることを特徴とする付記１に記載のプログラム。
【００８５】
（付記４）前記固有ベクトルを求めるステップにおいて、前記ハウスホルダー変換のための行列は、各プロセッサで処理できるような小さな係数行列となる上三角行列を各プロセッサで同時に生成されることを特徴とする付記１に記載のプログラム。
【００８６】
（付記５）前記固有ベクトルを求めるステップにおいて、前記固有ベクトルは、格納されている２次元配列の２次元目をプロセッサの数に対応して均等に分割し、おのおのの分割された領域を各プロセッサに割り当てることによって、演算されることを特徴とする付記１に記載のプログラム。
【００８７】
（付記６）共有メモリ型スカラ並列計算機用固有値問題の並列処理方法であって、
演算すべき実対称行列あるいはエルミート行列をブロック化し、ブロック化した部分をメモリの作業域にコピーして、ブロック間の積を利用して３重対角化するステップと、
３重対角化された行列から固有値と固有ベクトルを求めるステップと、
該３重対角化された行列から求められた固有ベクトルを、所定のブロック幅の行列演算の並列演算となるように変形されたハウスホルダー変換によって変換し、元の行列の固有ベクトルを求めるステップと、
を備えることを特徴とする並列処理方法。
【００８８】
（付記７）前記３重対角化ステップにおいて、ブロック化された部分の更新を再帰的プログラムで行うことを特徴とする付記６に記載の並列処理方法。
（付記８）前記３重対角化ステップにおいて、ブロック化された部分は、キャッシュメモリのページを跨らないように小さなブロックに分割して格納し、各プロセッサで並列に演算を行うことによって処理されることを特徴とする付記６に記載の並列処理方法。
【００８９】
（付記９）前記固有ベクトルを求めるステップにおいて、前記ハウスホルダー変換のための行列は、各プロセッサで処理できるような小さな係数行列となる上三角行列を各プロセッサで同時に生成されることを特徴とする付記６に記載の並列処理方法。
【００９０】
（付記１０）前記固有ベクトルを求めるステップにおいて、前記固有ベクトルは、格納されている２次元配列の２次元目をプロセッサの数に対応して均等に分割し、おのおのの分割された領域を各プロセッサに割り当てることによって、演算されることを特徴とする付記６に記載の並列処理方法。
【００９１】
【発明の効果】
高性能かつスケーラビリティのある固有値・固有ベクトルの解法を実現できる。
【図面の簡単な説明】
【図１】本発明の実施形態が前提とする共有メモリ型スカラ計算機のハードウェア構成を示した図である。
【図２】本発明の実施形態のアルゴリズムを説明する図（その１）である。
【図３】本発明の実施形態のアルゴリズムを説明する図（その２）である。
【図４】本発明の実施形態のアルゴリズムを説明する図（その３）である。
【図５】本発明の実施形態のアルゴリズムを説明する図（その４）である。
【図６】本発明の実施形態のアルゴリズムを説明する図（その５）である。
【図７】本発明の実施形態のアルゴリズムを説明する図（その６）である。
【図８】本発明の実施形態のアルゴリズムを説明する図（その７）である。
【図９】本発明の実施形態のアルゴリズムを説明する図（その８）である。
【図１０】本発明の実施形態のアルゴリズムを説明する図（その９）である。
【図１１】本発明の実施形態のアルゴリズムを説明する図（その１０）である。
【図１２】本発明の実施形態に従ったルーチンの疑似コード（その１）である。
【図１３】本発明の実施形態に従ったルーチンの疑似コード（その２）である。
【図１４】本発明の実施形態に従ったルーチンの疑似コード（その３）である。
【図１５】本発明の実施形態に従ったルーチンの疑似コード（その４）である。
【図１６】本発明の実施形態に従ったルーチンの疑似コード（その５）である。
【図１７】本発明の実施形態に従ったルーチンの疑似コード（その６）である。
【図１８】本発明の実施形態に従ったルーチンの疑似コード（その７）である。
【図１９】疑似コードの処理をフローチャートで示した図（その１）である。
【図２０】疑似コードの処理をフローチャートで示した図（その２）である。
【図２１】疑似コードの処理をフローチャートで示した図（その３）である。
【図２２】疑似コードの処理をフローチャートで示した図（その４）である。
【図２３】疑似コードの処理をフローチャートで示した図（その５）である。
【図２４】疑似コードの処理をフローチャートで示した図（その６）である。
【図２５】疑似コードの処理をフローチャートで示した図（その７）である。
【図２６】疑似コードの処理をフローチャートで示した図（その８）である。
【図２７】疑似コードの処理をフローチャートで示した図（その９）である。
【図２８】疑似コードの処理をフローチャートで示した図（その１０）である。
【図２９】疑似コードの処理をフローチャートで示した図（その１１）である。
【符号の説明】
１０−１〜１０−ｎ　　　プロセッサ
１１−１〜１１−ｎ　　　メモリモジュール
１２　　　相互結合網
１３−１〜１３−ｎ　　　２次キャッシュメモリ
【図１】
2004005528.tif 000008
【図２】
2004005528.tif 000009
【図３】
2004005528.tif 000010
【図４】
2004005528.tif 000011
【図５】
2004005528.tif 000012
【図６】
2004005528.tif 000013
【図７】
2004005528.tif 000014
【図８】
2004005528.tif 000015
【図９】
2004005528.tif 000016
【図１０】
2004005528.tif 000017
【図１１】
2004005528.tif 000018
【図１２】
2004005528.tif 000019
【図１３】
2004005528.tif 000020
【図１４】
2004005528.tif 000021
【図１５】
2004005528.tif 000022
【図１６】
2004005528.tif 000023
【図１７】
2004005528.tif 000024
【図１８】
2004005528.tif 000025
【図１９】
2004005528.tif 000026
【図２０】
2004005528.tif 000027
【図２１】
2004005528.tif 000028
【図２２】
2004005528.tif 000029
【図２３】
2004005528.tif 000030
【図２４】
2004005528.tif 000031
【図２５】
2004005528.tif 000032
【図２６】
2004005528.tif 000033
【図２７】
2004005528.tif 000034
【図２８】
2004005528.tif 000035
【図２９】
2004005528.tif 000036