JP4037303B2 - 共有メモリ型スカラ並列計算機用固有値問題の並列処理方法 - Google Patents

Description

【０００１】
【発明の属する技術分野】
本発明は、共有メモリ型スカラ並列計算機における行列の演算に関する。
【０００２】
【従来の技術】
ベクトル並列計算機用実対称行列（行列要素が実数で、転置しても行列が不変に保たれる行列）及びエルミート行列（行列要素が複素数で、共役・転置をとっても不変に保たれる行列）の固有値問題（Ａを行列、λを定数、Ｉを単位行列とした場合、det｜Ａ−λＩ｜＝０となるλとその固有ベクトルを求めること）は、３重対角化（対角成分とその両隣の成分のみを持つ行列にすること）をまず行う。この３重対角行列の固有値問題をマルチセクション法によって解き、固有値を求めてから逆反復法で固有ベクトルを求める。その後、固有ベクトルにハウスホルダー変換を行い、元の固有値問題の固有ベクトルを求めていた。
【０００３】
【発明が解決しようとする課題】
ベクトル計算機においては、メモリアクセス系の動作速度が速いことを前提として、固有値問題の演算を行っていたが、共有メモリ型スカラ計算機の場合、演算する行列が大きくなるほどに、共有メモリにアクセスする回数が多くなり、アクセス速度の遅い共有メモリへのアクセスによって、計算機の性能が大きく損なわれてしまうという問題がある。そこで、共有メモリ型スカラ計算機の各プロセッサに設けられる、アクセス速度の速いキャッシュメモリを有効に使い、行列計算を行う必要がある。つまり、行列の各行あるいは列毎に演算していると、共有メモリへのアクセスが多くなってしまうので、行列をブロック化して、各プロセッサにキャッシュメモリに格納されたデータを最大限処理した後、共有メモリにアクセスするようにして、共有メモリへのアクセス数を減少する、各プロセッサに局所化したアルゴリズムが必要となる。
【０００４】
すなわち、共有メモリ型の並列計算機では、ベクトル計算機のような強力なメモリアクセス能力がないため、メモリアクセスに対して演算量を増やすアルゴリズム上の工夫が必要になる。
【０００５】
本発明の課題は、共有メモリ型スカラ計算機において、高速に固有値問題の演算が行える演算の並列処理方法を提供することである。
【０００６】
【課題を解決するための手段】
本発明の並列処理方法は、共有メモリ型スカラ並列計算機用固有値問題の並列処理方法を実現させるプログラムであって、演算すべき実対称行列あるいはエルミート行列をブロック化し、ブロック化した部分をメモリの作業域にコピーして、ブロック間の積を利用して３重対角化するステップと、３重対角化された行列から固有値と固有ベクトルを求めるステップと、該３重対角化された行列から求められた固有ベクトルを、所定のブロック幅の行列演算の並列演算となるように変形されたハウスホルダー変換によって変換し、元の行列の固有ベクトルを求めるステップとを備えることを特徴とする。
【０００７】
本発明によれば、固有値問題を共有メモリ型スカラ計算機の各プロセッサにできるだけ局所化した演算で解くことができるので、各プロセッサが、共有メモリに頻繁にアクセスするために生じる遅延時間が最小限に抑えられ、並列計算の威力を有効に発揮した演算が行える。
【０００８】
【発明の実施の形態】
本発明の実施形態においては、上記固有値問題の演算における３重対角化部分のアルゴリズムでは、ブロック化されたアルゴリズムを採用し、ブロック部分を求める部分のアルゴリズムを再帰的に行うようにし、更新部分での演算密度を高める。また、メモリのアクセスがキャッシュメモリのページ間をとびとびにアクセスすることを避けるように行列ベクトル積部分を対称性を利用して連続アクセスが使える方式にする。キャッシュメモリのページを跨ってデータを読み出す場合、一度に読み込めず、キャッシュメモリに２度アクセスする必要が出てくるので、計算機のパフォーマンスを落とすことにある。従って、上記のように、ページをデータがまたがないようにする。
【０００９】
また、３重対角化された行列の固有ベクトルにハウスホルダー変換を行い、元の行列の固有ベクトルを求める部分は、ハウスホルダー変換を８０回程度ごとに束ねて、３つの行列要素からなる演算を行うことで演算密度を高める。
【００１０】
なお、本発明の実施形態においては、３重対角化された行列から固有値を求める方法及び３重対角化された行列の固有ベクトルを求める方法は、従来の方法によるとする。
【００１１】
図１は、本発明の実施形態が前提とする共有メモリ型スカラ計算機のハードウェア構成を示した図である。
プロセッサ１０−１〜１０−ｎは、１次キャッシュメモリを持っており、この１次キャッシュメモリはプロセッサの中に組み込まれていることもある。また、各プロセッサ１０−１〜１０−ｎには、２次キャッシュメモリ１３−１〜１３−ｎが設けられ、２次キャッシュメモリ１３−１〜１３−ｎが相互結合網１２に結合されている。また、相互結合網１２には、共有メモリであるメモリモジュール１１−１〜１１−ｎが設けられ、プロセッサ１０−１〜１０−ｎは、演算に必要なデータをここから読み出し、相互結合網１２を介して、２次キャッシュメモリ１３−１〜１３−ｎあるいは、１次キャッシュメモリに記憶させて、演算を行う。
【００１２】
この場合、メモリモジュール１１−１〜１１−ｎから２次キャッシュメモリ１３−１〜１３−ｎあるいは１次キャッシュメモリにデータを読み込んだり、２次キャッシュメモリ１３−１〜１３−ｎ、あるいは、１次キャッシュメモリから演算後のデータをメモリモジュール１１−１〜１１−ｎに書き込むのはプロセッサ１０−１〜１０−ｎの演算速度に比べて非常に遅い。従って、このような書き込み、読み出しが頻繁に発生すると計算機全体の性能を劣化させてしまう。
【００１３】
従って、計算機全体の性能を高く維持するためには、メモリモジュール１１−１〜１１−ｎへのアクセスをできるだけ減らし、２次キャッシュメモリ１３−１〜１３−ｎ、１次キャッシュメモリ、及びプロセッサ１０−１〜１０−ｎからなるローカルな系においてできるだけ多くの演算をしてしまうようなアルゴリズムが必要となる。
＜固有値、固有ベクトルを求める方式＞
１．３重対角化部分
１）３重対角化
ａ）ブロック化された３重対角化の数学的アルゴリズム
行列をブロック幅ごとに３重対角化する。行列をブロックに分割してブロック単位で下記のアルゴリズムで３重対角化を行う。
【００１４】
図２〜図１１は、本発明の実施形態のアルゴリズムを説明する図である。
図２は、ｍ番目のブロックを処理しているときの図である。ブロックとは図２で点線で示した列及び行方向の矩形である。
【００１５】
最後のブロックは、そのブロック幅−２をブロック幅とする、右下に２×２行列を残した左詰めのブロックに関してアルゴリズムを適用して全ての処理を終わる。
【００１６】
do i=1,blks
step1:A_nのn+1行ベクトルよりハウスホルダーベクトルuを作る。
step2:v_i＝A_n+iu、w_i＝v_i-u（u^tv）/2を計算する。
step3:U_i＝（U_i-1、u_i）、W_i＝（W_i-1、w_i）と更新する。（ここで、（Ｕ_i-1、ｕ_ｉ）は、行列Ｕ_i-1という１列を加えて行列Ｕ_iとすることにより、行列を１列づつ拡大する。）
step4: if(i<blks) then
A_nのn+i+1列目を更新する。
【００１７】

ｂ）ハウスホルダー変換を構成する情報の格納
固有ベクトルを計算するときに３重対角化で使ったハウスホルダー変換が必要になる。このためU_nおよびαをハウスホルダー変換を構成するベクトルの位置に格納する。αは対応する対角要素の位置に格納する。
ｃ）U_iを効率的に求める方法
ブロック部分の３重対角化を行うためにハウスホルダー変換のための次のベクトルを更新する必要がある。これらの計算をできるだけローカルに行うために、ブロック幅部分を作業域にコピーして３重対角化を行い、元の領域に格納する。次の列ベクトルの更新をそのたび行わず、行列積の形にして演算密度を上げて計算を行う。このため、ブロック部分の３重対角化を再帰的なプログラミングで行う。
【００１８】
recursive subroutine trid（幅、ブロック領域ポインタ）
if(幅＜１０）then
ｃ 当該幅のブロック部分の３重対角化を行う。
【００１９】
ハウスホルダー変換で必要なベクトルｕ、及び行列ベクトル積からｖ_iおよびｗ_iを作成する。
Ｕ、Ｗにｕ_i、ｗ_iを結合する。
【００２０】
else
ｃ ブロック幅を半分に分割
ｃ 前半ブロックに関して３重対角化を行う。
【００２１】
call trid（前半の幅、前半の領域）
ｃ ブロックを分割して、分割線から後半部分を更新
B=B-UW^t-WU^tで更新する。
ｃ その後、後半部分を３重対角化する。
【００２２】
call trid（後半の幅、後半の領域）
return
end
図３に示されるように、作業領域Ｕにブロック部分をコピーして、ブロック部分を再帰的プログラムで３重対角化する。再帰的プログラムのため、図３の前半部分は前半の更新処理のために当該再帰的プログラムが呼び出されたときに３重対角化される。後半部分は、前半部分によって更新され、その後、３重対角化される。
【００２３】
図４に示されるように、再帰的プログラムが深さ２まで呼ばれると、最初の前半処理で図４（ａ）の斜線部分をＢとして更新し、次に、図４（ｃ）の斜線部分、最後に図４（ｅ）が更新される。更新時の並列化は、更新部分のブロック行列を行ベクトル方向に均等に分割して、各部分の更新を複数のプロセッサで並列に行う。
図４（ａ）の演算後、図４（ｂ）の演算が行われ、図４（ｃ）の演算後、図４（ｄ）の演算が行われ、図４（ｅ）の演算後、図４（ｆ）の演算が行われる。
【００２４】
図５に示されるように、Ｕの斜線部分を更新するとき、ｕの横線及びＷの縦線部分を参照する。このことで、演算密度を上げることができる。すなわち、
A_n+k=A_n-U_kW_k ^t-W_kU_k ^t .......(**)
を利用して、V_nを次のように計算できる。このときU、Wの参照パターンは
v_k=A_nu_k-U_k-1W_k-1 ^tu_k-W_k-1U_k-1 ^tu_k .....(***)
図５（ａ）と（ｂ）、（ｃ）と（ｄ）、（ｅ）と（ｆ）のＵの更新の後、更新した部分の３重対角化のためにv_kを計算するときには、Ｕ、Wを参照して、行列ベクトル積で計算する。参照だけであり、Ｕの更新と参照が共通する部分もあり、効率的に参照できる。A_nを毎回更新せずに、このようにU、Wを使って必要な部分を更新し、（＊＊）を使って、全体の更新の演算密度を上げて性能向上を行う。（＊＊＊）は余分な演算ではあるが、ブロック幅を小さくしておく限り無視できる。
【００２５】
例えば、４個のＣＰＵで並列処理を行う場合、（＊＊＊）のW_k-1 ^tu_k、U_k-1 ^tu_kの計算では、縦線での分割で並列計算を行い、その結果との積に関しては、破線の方向に分割して並列計算する。
v_i=A_nu_iの計算の並列化
図６に示されるように、A_nの対称性、つまり、A_n＝A_n ^tを利用して各プロセッサで２次元目を分割して、v_iの部分をA_n（*、ns：ne）^tu_iを各プロセッサで計算する。
２）共有メモリ型スカラ並列計算機での並列化
ａ）U、W用の格納領域は、共有メモリに確保する。また、３重対角化を行うブロック領域は別に確保した作業領域にコピーして、その領域で３重対角化を行う。
【００２６】
上記の再帰的プログラムの中で、並列化する部分は以下の部分
（１）ハウスホルダー変換を構築するために必要なu_iを求めるためにstep4の
A_n(*, n+i+1)=A_n(*, n+i+1)-U_iW_i(n+i+1,*)^t-W_iU_i(n+i+1,*)^t
を使って、必要なベクトルを求める部分
（２）step2でv_iを計算する部分。この部分は（＊＊）の以下の式にu_iを作用させて計算する。
A_n+k=A_n-U_kW_k ^t-W_kU_k ^t
この計算ではA_nとu_iとの積及びU_kW_k ^t-W_kU_k ^tとu_iとの積を並列化する。
【００２７】
ブロック部分は、作業域にコピーされておりA_nの必要部分は更新されないようにする。また、A_nの対称性を利用して列ベクトル方向に長い行列に分割して並列計算する。
（３）再帰的プログラムで、ブロック領域の更新は、
A_n+k=A_n-U_kW_k ^t-W_kU_k ^t
を利用して更新する。そのことで（１）の計算を減らす。
３）step5の更新部分
更新は対称性を利用して対角成分より下半分を計算する。並列化は負荷がバランスするようにＣＰＵ数を＃ＣＰＵとしたとき、更新すべき部分行列の格納されているサブアレイの２次元目を２×＃ＣＰＵ個に均等に分割して、１から２×＃ＣＰＵと番号を振る。１から＃ＣＰＵのＣＰＵに関してｉ番目のプロセッサでｉ番目の部分と２×＃ＣＰＵ＋１−ｉ番目の部分の更新を並列に行う。
【００２８】
このあと上半分にコピーを行う。これも同様に分割して負荷をバランスさせる。このとき、対角ブロック部分以外は、キャッシュメモリのページをまたがらないようにある程度小さなブロックに分割してコピーを行う。
A_n+k=A_n-U_kW_k ^t-W_kU_k ^tで更新する部分は下三角行列部分で、図７のように列ブロックにＣＰＵ数×２に分割して、両端からペアを作り、これを各ＣＰＵで更新する。図７は、４個のＣＰＵが設けられている場合を記述している。
【００２９】
下三角部分が更新できたら、上三角部分に１〜８を同じペアで転置してｕ１からｕ８にコピーする。
このとき、小さな正方ブロックに内部を分割してキャッシュを利用して転置する。更新と同様にＣＰＵで並列処理する。
キャッシュ上での転置で性能向上する部分の説明
図８に示されるように、１から番号順に転置変換していく。各正方形領域に関して、１の場合は下三角部分をメモリの連続領域にコピーして、それを行方向にアクセスして（転置して）１の上三角部分に格納する。２〜８は正方形にコピーして、行方向に参照して列方向（上下方向）に格納領域に格納する。
２固有ベクトルの計算
ａ）基本的なアルゴリズム
ベクトルu_nを格納して、後で（1-2*uu^t/(u^tu））を構成して、ベクトルにかける。
【００３０】
３重対角化を行うと元の固有値問題は以下のようになる。
Q_n-2......Q₂Q₁AQ₁ ^tQ₂ ^t........Q_n-2 ^tQ_n-2......Q₂Q₁x=λQ_n-2......Q₂Q₁x
３重対角化した固有値問題を解いて求めた固有ベクトルｙからx=Q₁ ^tQ₂ ^t.....Q_n-3 ^tQ_n-2 ^tyを計算することで変換することができる。
ｂ）本発明の実施形態のブロックアルゴリズム及び固有ベクトルの変換の並列計算
多数のまたは全ての固有ベクトルを求める場合は各ＣＰＵに均等に３重対角行列の固有ベクトルを対応付けて各ＣＰＵで並列に上記の変換を施す。このとき変換行列を８０本程度づつまとめて変換を行う。
【００３１】
各々の変換行列Q_i ^tは1+α_iu_iu_i ^tと書ける。これらの積を展開すると以下の形になる。
b_i,jは最左端及び最右端のu_iu_jを除いたスカラ係数をまとめたものである。
【００３２】
【数１】

【００３３】
の形になる。b_i,jは上三角行列になる。これを求めて１＋ＵＢＵ^tの形になる。このような形に変形することで、計算密度を向上させることで高速化を図る。図９は、行列Ｂの模式図である。
【００３４】
上記の方法は３ステップであるが、メモリアクセスの関係から行列はＵ及びＢの２つとなり、Ｂは比較的小さな領域にすることができるため効率がよい。
ｍ−１番目までb_i,jができた後に、（１＋α_mＵ_mＵ_m ^t）をかけると、
【００３５】
【数２】

【００３６】
最後の項の和をｉ、ｊを入れ替えて
【００３７】
【数３】

【００３８】
と変更することができる。
最内の括弧の中がb_m,j（j=m+1、・・・、n+k）と考えられる。b_m,mはα_mである。
ブロック幅の正方な作業配列Ｗ２を用意して、はじめにα₁U_iU_j ^tをｗ２（ｉ、ｊ）の上三角部分に格納する。対角要素にはα_iを格納する。
【００３９】
次に、以上の議論は最右下の２×２上三角行列から順番に１行ずつ上に付け加える要領で計算することができる。
行の最も右の要素から一つずつ計算していけば、Ｂは上三角行列であるため更新された部分の参照が無いために同じ領域で計算ができる。このように、極めて小さな領域のみを使って３つの行列積の真ん中に位置する係数行列を計算することができる。
【００４０】
図１０は、上記固有値を求める方法を模式的に示した図である。
ブロック幅をnbsとする。
まず、Ｂの上半分に内積α_jｕ_i・ｕ_jを計算して格納する。
【００４１】
対角要素にはα_iを格納する。
次に、以下のように計算する。
do i1=nbs-2,1,-1
do i2=nbs,i1+1,-1
sum=w2(i1,i2)
do i3=i2-1, i1+1, -1
sum=sum+w2(i1,i3)*w2(i3, i2)
enddo
w2(i1,i2)=sum
enddo
enddo
do i2=nbs,1,-1
do i1=i2-1,1,-1
w2(i1,i2)=w2(i1,i2)*w2(i2,i2)
enddo
enddo
図１１は、上記で求められた固有ベクトルを元の行列の固有ベクトルに変換する様子を模式的に示した図である。
【００４２】
配列Ａに格納されているハウスホルダーベクトルで固有ベクトルを変換する。ブロック化されていて、図１１の斜線の部分をＥＶの斜線部分にかけて結果をＷに格納する。また、Ａブロック行列からＷ２を構成する。Ｗ２とＷをかけたものにＡのブロック部分をかけたものからＥＶの斜線部分を更新する。
３エルミート行列の固有値・固有ベクトル
エルミート行列の固有値・固有ベクトルを求めるアルゴリズムは実対称行列の３重対角化での転置を転置＋複素共役に置き換える（ｔ→Ｈ）。また、ハウスホルダーベクトルの作り方は最初の要素のスカラ倍になるようにベクトルの大きさを加える。
【００４３】
できた３重対角行列はエルミート行列であり、これを絶対値１の対角行列でスケーリングする。
対角行列は以下のように作る。
ｄ_i＝１．０、ｄ_i+1＝ｈ_i+1,i/｜ｈ_i+1｜＊ｄ_i
図１２〜図１８は、本発明の実施形態に従ったルーチンの疑似コードである。
【００４４】
図１２は、実対称行列を３重対角化するサブルーチンである。
配列ａに実対称行列の下三角部分を格納する。ｄａｉｇ、ｓｄｉａｇに３重対角行列、副対角部分を格納する。出力としてａの下三角部分に変換で必要な情報を格納する。
【００４５】
Ｕは３重化を行うブロックを格納する。ＶはＷを格納する領域である。
ｎｂは、ブロックの数であり、ｎｂａｓｅはブロックの始まり位置を示す。
サブルーチンｃｏｐｙを実行した後、３重対角化をするブロックをｕ（ｎｂａｓｅ＋１：ｎ、１：ｉｂｌｋ）に格納し、ｂｌｋｔｒｉｄルーチンを呼び出してＬＵ分解を行う。そして、処理の終わったｕ（ｎｂａｓｅ＋１：ｎ、１：ｉｂｌｋ）を元の行列ａに書き戻す。これらの下の処理は、最後の残りのブロックの３重対角化をサブルーチンｂｌｋｔｒｉｄを用いて行うものである。
【００４６】
図１３は、３重対角化を行うサブルーチンの疑似コードである。
このサブルーチンは、ブロック行列を３重対角化するルーチンで、再帰的に呼び出される。ｎｂａｓｅはブロックの位置を示すオフセット、ｉｓｔａｒｔは再帰呼び出しで対象となる縮小されたサブブロックのブロック内でのオフセットで最初１で再帰的に呼び出されるとき対象サブブロックの位置を示す。ｎｗｉｄｔｈはサブブロックの大きさである。
【００４７】
ｎｗｉｄｔｈが１０より小さい場合には、サブルーチンｂｔｕｎｉｔを呼び出し、そうでない場合には、ｉｓｔａｒｔをｉｓｔａｒｔ２に格納し、ｎｗｉｄｔｈの２分の１をｎｗｉｄｔｈ２に格納し、サブルーチンｂｌｋｔｒｉｄによって３重対角化し、その後、バリア同期を取る。
【００４８】
更に、ｉｓｔａｒｔとｎｗｉｄｔｈ／２を加えたものをｉｓｔａｒｔ３に格納し、ｎｗｉｄｔｈ−ｎｗｉｄｔｈ／２をｎｗｉｄｔｈ３に格納し、ブロックの始まりと終わりをそれぞれ示すｉｓ２、ｉｓ３、ｉｅ２、ｉｅ３、ｉｓ、ｉｅに値を設定し、ｌｅｎ、ｉｐｔｒを設定し、同図の式に従って計算した後、ｕ（ｉｓ：ｉｅ、ｉｓ３：ｉｅ３）に格納して、バリア同期を取り、３重対角化サブルーチンｂｌｋｔｒｉｄを呼び出して処理をしてサブルーチン処理を終わる。
【００４９】
図１４は、３重対角化サブルーチンの内部ルーチンの疑似コードである。
３重対角化の内部サブルーチンｂｔｕｎｉｔにおいては、必要な情報を格納した後、ブロックの始まりｉｐｔｒ２と、幅ｌｅｎ、始まり位置ｉｓ、終わり位置ｉｅを決定し、バリア同期を取る。そして、ｕ（ｉｓ：ｉｅ、ｉ）^ｔ＊う（ｉｓ：ｉｅ、ｉ）をｔｍｐに格納し、バリア同期を取る。その後、各値を計算し、それぞれの配列に格納する。ここでｓｕｍは和、ｓｑｒｔは平方根をとることを意味する。そして、最後にバリア同期を取る。
【００５０】
次に、ｖ（ｉｓ：ｉｅ、ｉ）を求めて、バリア同期を取り、ｌｅｎ２、ｉｓｘ、ｉｅｘ、ｕ、ｖを更新して、バリア同期を取る。更に、ｖ（ｉｓ：ｉｅ、ｉ）を更新し、バリア同期を取る。更に、ｖ（ｉｓ：ｉｅ、ｉ）^t＊ｕ（ｉｓ：ｉｅ、ｉ）を計算し、ｔｍｐに格納し、バリア同期を取る。
【００５１】
そして、ｂｅｔａに値を設定し、バリア同期を取り、ｖをｂｅｔａを使った演算で更新し、バリア同期を取る。
そして、ｉ＜ｉｂｌｋかつｐｔｒ２＜ｎ−２である時には、ｕ（ｉｓ：ｉｅ、ｉ＋１）を更新し、他の場合には、別の式でｕ（ｉｓ：ｉｅ、ｉ；１：ｉ＋２）を更新し、処理を終了する。なお、このサブルーチンの実行後、スレッドの割当てを解除する。
【００５２】
図１５は、ｕ、ｖから行列の下半分を更新するルーチン及び、対角行列部分の更新をするルーチン及びコピーのルーチンの疑似コードである。
ここで、ｎｂａｓｅはブロックの位置を示すオフセットであり、ｎｗｉｄｔｈは、ブロック幅である。
【００５３】
このサブルーチンｕｐｄａｔｅでは、配列ａ、ｕ、ｖを確保した後、バリア同期を取る。そして、ｂｌｋ、ｎｂａｓｅ２、ｌｅｎ、ｉｓ１、ｉｅ１、ｎｂａｓｅ３、ｉｓｒ、ｉｅｒを設定した後、ａ（ｉｅ１：ｎ、ｉｓ１：ｉｅ１）、及び、ａ（ｉｅｒ＋１：ｎ、ｉｓｒ：ｉｅｒ）をそれぞれ更新し、サブルーチンｔｒｕｐｄａｔｅを２回呼び出し、バリア同期を取って処理を返す。
【００５４】
サブルーチンｔｒｕｐｄａｔｅでは、対角ブロック更新用ブロック幅をｂｌｋ２に格納し、ａ、ｕ、ｖの配列を確保する。
ｉｓ２、ｉｅ２を設定し、Ａ（ｉｓ２：ｉｅ、ｉｓ２：ｉｅ２）を更新し、処理を返す。
【００５５】
サブルーチンｃｏｐｙでは、ｌｅｎ、ｉｓ１、ｌｅｎ１、ｎｂａｓｅ、ｉｓｒ、ｌｅｎｒを設定して、ｂａｎｄｃｐを２回実行し、処理を返す。
図１６は、更新した下三角部分を上三角部分にコピーするルーチンの疑似コードである。
【００５６】
サブルーチンｂａｎｄｃｐでは、ｎｂ、ｗ、ｎｎ、ｌｏｏｐｘを設定し、ｄｏループの中で、ＴＲＬ（ｗ（１：ｎｎｘ、１：ｎｎｘ））にＴＲＬ（ａ（ｉｓ２：ｉｓ２＋ｎｎｘ−１、ｉｓ２：ｉｓ２＋ｎｎｘ））を格納し、ＴＲＵ（ａ（ｉｓ２：ｉｓ２＋ｎｎｘ−１、ｉｓ：ｉｓ＋ｎｎｘ））にＴＲＬ（ｗ（ｉ：ｎｎｘ、１：ｎｎｘ））^tを格納する。ここで、ＴＲＬは下三角部分、ＴＲＵは上三角部分を示す。
そして、ｗ（１：ｎｎｘ、１：ｎｎｘ）とａ（ｉｓ２：ｉｓ２＋ｎｎｘ、ｉｓ３：ｉｓ３＋ｎｎｘ−１）を更新し、次に、ｗ（１：ｎｙ、１：ｎｘ）とａ（ｉｓ２：ｉｓ２＋ｎｎｘ、ｉｓ３：ｎ）を更新する。
【００５７】
そして、ｄｏループを抜けると、処理を返す。
図１７は、３重対角行列の固有ベクトルを元の行列の固有ベクトルに変換するルーチンの疑似コードである。
【００５８】
ここで、３重対角行列の固有ベクトルは、ｅｖ（１：ｎ、１：ｎｅｖ）に格納されている。また、ａは、３重対角化の出力であり、変換に必要な情報を下三角部分に格納している。
【００５９】
サブルーチンｃｏｎｖｅｖは、配列としてａとｅｖを取る。
スレッドを構成し、処理を並列処理で行う。
スレッド間でバリア同期をとり、ｌｅｎ、ｉｓ、ｉｅ、ｎｅｖｔｈｒｄを設定し、ｃｏｎｖｅｖｔｈｒｄルーチンを呼び出し、リターン後バリア同期を取って処理を終了する。
【００６０】
図１８は、固有ベクトルを変換するルーチンの疑似コードである。
サブルーチンｃｏｎｖｅｖｔｈｒｄでは、ブロック幅をｂｌｋに格納し、配列としてａ、ｅｖ、ｗ、ｗ２を取る。
【００６１】
まず、ｗｉｄｔｈが０より小さければ、処理をせず元のルーチンに戻る。ｎｕｍｂｌｋとｎｆｂｓを設定し、ａｌｐｈａに３重対角化時に対角要素に格納された値の符号を逆にしたもの（−ａ（ｉ、ｉ））を入れる。また、ｘ（１：ｉｗｉｄｔｈ）にｅｖ（ｉ＋１：ｎ、１：ｉｗｉｄｔｈ）^t＊ａ（ｉ＋１：ｎ、ｉ）を入力し、これと、ａｌｐｈａとａを用いてｅｖを更新する。更に、次のｄｏ文で、ｉｓ、ｉｅを設定し、ｗ（１：ｂｌｋ、ｉｗｉｄｔｈ）をａ（ｉｓ＋１：ｎ、ｉｓ：ｉｅ）^t＊ｅｖ（ｉｓ＋１：ｎ、１：ｉｗｉｄｔｈ）で上書きし、同じくｗ（１：ｂｌｋ、１：ｉｗｉｄｔｈ）をＴＲＬ（ａ（ｉｅ＋１：ｉｓ、ｉｓ：ｉｅ））^t＊ｅｖ（ｉｅ＋１：ｉｓ、１：ｉｗｉｄｔｈ）で更新する。ここで、ＴＲＬは、下三角行列部分である。
【００６２】
また、ｗ２の対角要素ベクトルＤＩＡＧ（ｗ２）に、ａ（ｉｓ：ｉｅ、ｉｓ：ｉｅ）の対角要素ベクトルを格納する。
次のｄｏ文でｗ２（ｉ１、ｉ２）をｗ２（ｉ１、ｉ２）＊（ａ（ｉｓ＋ｉ２：ｎ、ｉｓ＋ｉ２−１）^t＊ａ（ｉｓ＋ｉ２：ｎ、ｉｓ＋ｉ１−１））で更新する。更に次のｄｏ文で、ｗ２（ｉ１、ｉ２）をｗ２（ｉ１、ｉ２）＋ｗ２（ｉ１、ｉ１＋１：ｉ２−１）＊ｗ２（ｉ１＋１：ｉ２−１、ｉ２）で更新する。
【００６３】
更に、次のｄｏ文で、ｗ２（ｉ１、ｉ２）をｗ２（ｉ１、ｉ２）＊ｗ２（ｉ２、ｉ２）で更新し、ｗ（１：ｂｌｋ、１：ｉｗｉｄｔｈ）、ｅｖ（ｉｓ＋ｎ：ｎ、１：ｉｗｉｄｔｈ）、ｅｖ（ｉｅ＋１：ｉｓ、１：ｉｗｉｄｔｈ）を更新して、処理を元のルーチンに返す。
【００６４】
図１９〜図２９は、疑似コードの処理をフローチャートで示した図である。
図１９は、実対称行列を３重対角化するサブルーチンtridのフローである。ステップＳ１０においては、サブルーチンとしてshared配列Ａ（ｋ、ｎ）、diag（ｎ）およびsdiag（ｎ）を入力。diag、sdiagは計算した３重対角行列の対角要素、副対角要素を出力として返却する。作業域Ｕ（ｎ＋１、iblk）、ｖ（ｎ＋１、iblk）をルーチン内部で確保しshared属性で利用する。ステップＳ１１において、スレッド生成し、各スレッドでローカル域numthrに総スレッド数nothrdに各スレッドに割り振られたスレッド番号を設定し、各スレッドで以下を設定する。iblkにブロック幅を設定し、nb=(n-2+iblk-1)/iblk、nbase=0、i=1を設定する。ステップＳ１２においては、i>nb-1か否かを判断する。ステップＳ１２の判断がＹＥＳの場合には、ステップＳ１９に進む。ステップＳ１２における判断がＮＯの場合には、ステップＳ１３において、nbase=(i-1)×iblk、istart=1、nwidth=iblkを設定する。ステップＳ１４において、サブルーチンcopyを呼んで上三角部分に下三角部分をコピーする。ステップＳ１５において、作業領域Ｕにブロック三重対角化の対象領域をコピーする。すなわち、Ｕ（nbase+1:n、1:iblk）←Ａ（nbase+1:n、nbase+1:nbase+iblk）を実行する。ステップＳ１６において、サブルーチンblktridを呼んでＵにコピーした部分の三重対角化を行う。（istart=1、ブロック幅はiblkを受け渡す）。ステップＳ１７において、三重対角化されたものを配列Ａに戻す。すなわち、Ａ（nbase+1:n,nbase+1:nbase+iblk）←Ｕ（nbase+1:n、1:iblk）を行う。ステップＳ１８において、サブルーチンupdateを呼んで、Ａ（nbase+iblk:n、nbase+iblk:n）の下三角部分を更新し、ステップＳ１２に戻る。
【００６５】
ステップＳ１９においては、nbase=(nb-1)×iblk、istart=1、iblk2=n-nbaseと設定する。ステップＳ２０においては、作業領域Ｕにブロック三重対角化の対象領域をコピーする。すなわち、Ｕ（nbase+1:n,1:nwidth）←Ａ（nbase+1:n、nbase+1:n）を行う。ステップＳ２１において、サブルーチンblktridを呼んでＵにコピーした部分の三重対角化を行う（istart=1、ブロック幅はiblk2を受け渡す）。ステップＳ２２において、三重対角化されたものを配列Ａに戻す。すなわち、Ａ（nbase+1:n、nbase+1:n）←Ｕ（nbase+1:n、1:nwidth）を行う。ステップＳ２３においては、並列処理のために生成したスレッドを消して、サブルーチンを終了する。
【００６６】
図２０は、サブルーチンblktridのフローである。このサブルーチンは再帰的プログラムとなっている。
このサブルーチンは以下のステートメントによって呼び出される。
subroutine blktrid（A,k,n,dig,sdig,nbase,istart,nwidth、U、V、nothrd、numthrd）
ここで、nbaseはブロックの位置を示すオフセット。istartは再帰呼び出しで対称となる縮小されたブロックのブロック内でのオフセットで最初は１で再帰的に呼び出されるとき対象ブロックの位置を示す。nwidthはブロック幅を表す。ステップＳ２５においては、nwidth<10か否かを判断する。ステップＳ２５の判断がＮＯの場合には、ステップＳ２７に進む。ステップＳ２５の判断がＹＥＳの場合には、ステップＳ２６において、サブルーチンbtunitを呼んで３重対角化を行い、サブルーチンを抜ける。ステップＳ２７では、再帰呼び出しのために対象となる更新位置とブロック幅を変えて、呼び出すistart2=istart、nwidth2=nwidth/2を設定し受け渡す。縮小されたブロックの開始位置、ブロック幅を受け渡す。ステップＳ２８においては、サブルーチンblktridを再帰的に呼び出す。ステップＳ２９においては、スレッド間でバリア同期を取る。ステップＳ３０においては、更新で各スレッドが分担する始点（is2、is3）、終点（ie2、ie3）を計算する。すなわち、istart3=istart+nwidth/2、nwidth3=nwidth-nwidth/2、is2=istart2、ie2=istart+nwidth2-1、is3=istart3、ie3=istart3+nwidth3-1、iptr=nbase+istart3、len=(n-iptr+numthrd-1)/numthrd、is=iptr+(nothrd-1)×len+1、ie=min(n,iptr+nothrd×len）を計算する。ステップＳ３１においては、Ｕ（is:ie,is3:ie3）＝Ｕ（is:ie,is3:ie3）-Ｕ（is:ie,is2:ie2）×Ｗ（is3:ie3,is2:ie2）^t-Ｗ（is:ie,is2:ie2）×Ｕ（is3:ie3,is2:ie2）^tを計算し、ステップＳ３２において、スレッド間でバリア同期を取り、ステップＳ３３においてサブルーチンblktridを再帰的に呼び出して、サブルーチンを抜ける。
【００６７】
図２１及び図２２は、サブルーチンblktridの内部ルーチンであるサブルーチンbtunitのフローである。
ステップＳ３５において、tmp（numthrd）、sigma、alphaをshared属性で割り付ける。ステップＳ３６において、nbase+istart>n-2か否かを判断する。ステップＳ３６の判断がＹＥＳの場合には、サブルーチンを抜ける。ステップＳ３６の判断がＮＯの場合にはステップＳ３８に進む。ステップＳ３８では、i=istartとし、ステップＳ３９において、i<=istart-1+nwidthであるか否かを判断する。ステップＳ３９の判断がＮＯの場合には、サブルーチンを抜ける。ステップＳ３９の判断がＹＥＳの場合には、ステップＳ４０において、各スレッドで分担する始点is、終点ieを計算する。iptr2=nbase+i、len=(n-iptr2+numthrd-1)/numthrd、is=iptr2+(nothrd-1)×len+1、ie=min(n,iptr2+nothrd×len）を計算する。ステップＳ４１では、バリア同期を取り、ステップＳ４２では、tmp（nothrd）=Ｕ（is:ie,i）^t×Ｕ（is:ie,i）を計算し、ステップＳ４３で、バリア同期を取る。ステップＳ４４において、nothrd=1か否かを判断する。ステップＳ４４の判断がＮＯの場合には、ステップＳ４６に進む。ステップＳ４４の判断がＹＥＳの場合には、ステップＳ４５において、各スレッドで部分計算した和の平方根をとり、三重対角化のための計算を行う（ハウスホルダーベクトルの作成）。sigma=sqrt(sum(tmp(1:numthrd)))、ここで、ＳＵＭは和、ｓｑｒｔは平方根である。diag（iptr2）=u(iptr2,i)、sdiag(iptr2)=-sigma、Ｕ（nbase+i+1、i）＝Ｕ（nbase+i+1,i）+sign(u(nbase+i+1,i)×sigma、alpha=1.0/(sigma×u(nbase+i+1,i）、Ｕ（iptr2,i）=alphaを計算して、ステップＳ４６に進む。ステップＳ４６では、バリア同期を取り、ステップＳ４７では、iptr3=iptr2+1、ステップＳ４８では、Ｖ（is:ie,i）=Ａ（iptr3:n、iptr2+is:iptr2+ie）^tＵ（ptr3:n,i）を演算し、ステップＳ４９では、バリア同期を取る。
【００６８】
ステップＳ５０では、Ｖ（is:ie,i）=alpha×Ｖ（is:ie,i）-Ｖ（is:ie,1:i-1）×（Ｕ（iptr3:n,1:i-1）^t×Ｕ（iptr3:n、i））-Ｕ（is:ie,1:i-1）×（Ｖ（iptr3:n,1:i-1)^t×Ｕ（iptr3:n,i））を計算し、ステップＳ５１において、バリア同期を取る。ステップＳ５２においては、tmp（nothrd）=Ｖ（is:ie,i）^t×Ｕ（is:ie,i）を計算し、ステップＳ５３において、バリア同期を取る。ステップＳ５４においては、nothrd=1か否かを判断する。ステップＳ５４の判断がＮＯの場合には、ステップＳ５６に進み、ステップＳ５４の判断がＹＥＳの場合には、ステップＳ５５に進む。ステップＳ５５においては、beta=0.5×alpha×sum(tmp(1:numthrd))を計算する。ただし、sumはベクトルの和をとる記号である。ステップＳ５６では、バリア同期を取る。ステップＳ５７では、Ｖ（is:ie、i）=Ｖ(is:ie,i）-beta×Ｕ（is:ie,i）を計算し、ステップＳ５８において、バリア同期を取る。ステップＳ５９においては、ptr2<n-2か否かを判断し、判断がＹＥＳの場合には、ステップＳ６０において、Ｕ（is:ie,i+1）=Ｕ（is:ie,i+1）-Ｕ（is:ie,istart:i）×Ｖ（i+1,istart:1）^t-Ｖ（is:ie,istart:i）×Ｕ（n+1,istart:1）^tを計算し、ステップＳ３９に戻る。ステップＳ５９の判断がＮＯの場合には、ステップＳ６１において、Ｕ（is:ie,i+1:i+2）=Ｕ（is:ie,i+1:I+2）-Ｕ（is:ie,istart:i）×Ｖ（i+1:n,istart:i）^t-Ｖ（is:ie,istart:i）×Ｕ（n+1:n,istart:i）^tを計算して、サブルーチンを抜ける。
【００６９】
図２３は、サブルーチンupdateのフローである。
ステップＳ６５において、バリア同期を取り、ステップＳ６６において、各スレッドでペアを作り、更新を分担する始点、終点を決める。すなわち、nbase2=nbase+iblk、len=(n-nbase2+2×numthrd-1）/（2×numthrd）、is1=nbase2+(nothrd-1)len+1、ie1=min(n,nbase2+nothrd×len）、nbase3=nbase2+2×numthrd×len、isr=nbase3-nothrd×len+1、ier=min(n,isr+len-1）を計算する。ステップＳ６７では、Ａ（ie1+1:n,is1:ie1）=Ａ（ie1+1:n,is1+1:n,is1:ie1)-Ｗ（ie1+1:n,1:blk）×Ｕ（is1:ie1,1:blk）^t-Ｕ（ie1+1:n,1:blk）×Ｗ（is1:ie1,1:blk）^t、Ａ（ier+1:n,isr:ier）=Ａ（ier+1:n,isr:ier）-Ｗ（ier+1:n,1:blk）×Ｕ（isr:ier,1:blk）^t-Ｕ（ier+1:n,1:blk）×Ｗ（isr:ier,1:blk）^tを計算し、ステップＳ６８で、サブルーチンtrupdateを呼んで左半分の対角行列部分を更新する。ｉｓ１、ｉｅ１、Ａ、Ｗ、Ｕを渡す。ステップＳ６９において、サブルーチンtrupdateを呼んで右半分の対角行列部分を更新する。ｉｓｒ、ｉｅｒ、Ａ、Ｗ、Ｕを渡す。ステップＳ７０において、バリア同期を取って、サブルーチンを抜ける。
【００７０】
図２４は、サブルーチンtrupdate（対角行列部分の更新）のフローである。
更新の始点is、と終点ieを入力とする。これは、サブルーチンが呼び出される前に対角ブロックの下の長方形の部分の更新に使われた始点と終点である。
【００７１】
ステップＳ７５においてblk2に対角ブロック更新用ブロック幅を設定し、i=isを設定する。ステップＳ７６において、i>ie-1か否かを判断する。ステップＳ７６における判断がＹＥＳの場合には、サブルーチンを抜ける。ステップＳ７６における判断がＮＯの場合には、ステップＳ７７において、各スレッドでの更新の始点、終点を決める。すなわち、is2=i、ie2=min(i+blk2-1,ie-1）、Ａ（is2:ie-1,is2,ie2）=Ａ（is2:ie-1,is2,ie2）-Ｕ（is2:ie-1,1:blk）×Ｗ（is2:ie2,1:blk）^t-Ｗ（is2:ie-1,1:blk）×Ｕ（is2:ie2,1:blk）^tを計算し、ステップＳ７８において、i=i+blk2と設定し、ステップＳ７６に戻る。
【００７２】
図２５は、サブルーチンcopyのフローである。
ステップＳ８０において、コピーを各スレッドで対を作って並列実行するときの始点、と幅を計算する。すなわち、len=(n-nbase+2×numthrd-1）/（2×numthrd）、is1=nbase+（nothrd-1)×len+1、len1=max（0,min(n-is1+1,len））、nbase3=nbase+2×numthrd×len、isr=nbase3-nothrd×len+1、lenr=max（0,min(n-isr+1,len））を計算する。ステップＳ８１において、サブルーチンbandcpを呼び出す。ペアの左側の始点is1、幅len1できまる領域をコピーする。ステップＳ８２において、サブルーチンbandcpを呼び出す。ペアの右側の始点isr、幅lenrできまる領域をコピーする。
【００７３】
図２６は、サブルーチンbandcpのフローである。
このルーチンは、小さな作業域ＷＸを利用してキャッシュ上での転置を行いながらコピーを行う。作業域をＷＸ（nb、nb）とし、始点をis、幅をlenで受ける。
【００７４】
ステップＳ８５において、nn=min(nb,len）、loopx=（len+nn-1）/nn、j=1を演算する。ステップＳ８６において、j>loopxか否かを判断する。ステップＳ８６の判断がＹＥＳの時は、サブルーチンを終了する。ステップＳ８６の判断がＮＯの場合には、ステップＳ８７において、ＷＸにコピーする対角ブロックの大きさnnxとそのオフセットipを決める。ip=is+(j-1)×nn、n1=len-(j-1)×nn、nnx=min(nn,n1)、len2=n-ip-nnx+1、loopy=(len2+nn-1)/nn、TRL（WX（1:nnx,1:nnx））=TRL（Ａ（ip:ip+nnx-1,ip:ip+nnx-1））、TRU（Ａ（ip:ip+nnx-1,ip:ip+nnx-1））=TRL（WX(1:nnx,i:nnx））、i=1,is2=ip,is3=ip+nnxを計算する。ただし、TRUは上三角部分、TRLは下三角部分。
【００７５】
ステップＳ８８では、i>loopy-1か否かを判断する。ステップＳ８８の判断がＮＯの場合には、ステップＳ８９において、nn×nnxの部分を転置してコピーする。すなわち、ＷＸ（1:nn、1:nnx）＝Ａ（is3:is3+nn-1、is2:is2+nnx-1）、Ａ（is2:is2+nnx-1,is3:is3+nn-1）=WX（1,nn:1,nnx）^t、is3=is3+nnを計算して、ステップＳ８８に戻る。ステップＳ８８の判断がＹＥＳの場合には、ステップＳ９０で、最後の部分をコピーする。すなわち、nn=n-is3+1、ＷＸ（1:nn,1:nx）=Ａ（is3:n,is2:is2+nnx-1）、Ａ（is2:is2+nnx-1,is3:n）=ＷＸ（1:nn,1:nx）を演算して、ステップＳ８６に戻る。
【００７６】
図２７は、サブルーチンconvevのフローである。
このルーチンでは、求める固有ベクトルの数nev、ハウスホルダーベクトルはａの下半分に格納されている。３重対角行列の固有ベクトルがev(k,nev）に格納されている。
【００７７】
ステップＳ９５において、スレッドを生成する。各スレッドのローカル域のnumthr及びnothrdに総スレッド数とスレッド番号（１〜numthrd）を設定する。ステップＳ９６において、バリア同期を取る。ステップＳ９７において、各スレッドで分担計算する始点、終点を決める。すなわち、len=（nev+numthrd-1）/numthrd、is=(nothrd-1)×len+1、ie=min(nev,nothrd×len）、width=ie-is+1を計算する。ステップＳ９８において、サブルーチンconvevthrdを呼び出し、三重対角行列からもとの行列の固有ベクトルに変換する。各スレッドが担当する固有ベクトルの格納されている領域と固有ベクトルの数widthを渡す。ステップＳ９９において、バリア同期をとり、ステップＳ１００において、生成したスレッドを消滅させ、サブルーチンを抜ける。
【００７８】
図２８及び図２９は、サブルーチンconvevthrdのフローである。
このルーチンは、各スレッドごとに分担した三重対角行列の固有ベクトルを変換して、もとの行列の固有ベクトルにする。ハウスホルダ変換を復元するベクトルと係数は配列Ａに格納されている。
【００７９】
ステップＳ１１０において、ｂｌｋにブロック幅を設定する。ブロック幅は８０程度である。ステップＳ１１１では、iwidth<0か否かを判断する。ステップＳ１１１の判断がＹＥＳの場合には、サブルーチンを抜ける。ステップＳ１１１の判断がＮＯの場合には、ステップＳ１１２に進む。ステップＳ１１２では、以下のループで変換の最初のブロックは、（１＋αｕｕ^t）を順次計算して求める。まず、numblk=(n-2+blk-1)/blk、nfbs=n-2-blk×(numblk-1)を求める。ステップＳ１１３において、i<n-2-nfbs+1か否かを判断する。ステップＳ１１３の判断がＹＥＳの場合には、ステップＳ１１４に進み、alpha=-a(i,i)、x(1:iwidth)=a(i+1:n,i)^t×ev(i+1:n,1:width)、ev(i+1:n,1:width)=ev(i+1:n,1:width)+alpha×a(i+1:n,i)×x(1:iwidth)^tを演算し、ステップＳ１１３に戻る。ステップＳ１１３の判断がＮＯの場合には、ステップＳ１１５において、i=1に設定し、ステップＳ１１６において、i>numblk-1か否かを判断する。ステップＳ１１６の判断がＮＯの場合には、サブルーチンを抜ける。ステップＳ１１６の判断がＹＥＳの場合には、ステップＳ１１７において、ブロック形式の（１＋ＵＢＵ^t）のＵ^t×ＥＶをＵ^tの左端の上三角行列部分と右側の長方形の部分に分けて計算する。すなわち、is=n-2-(nfns+i×blk）+1、ie=ie+blk-1、W(1:blk,iwidth)=a(ie+1:n,is:ie)^t×ev(ie+1:n,1:iwidth)、W（1:blk-1,1:iwidth)=w(1:blk-1,1:iwidth）+TRL(a(is+1:ie,is:ie-1))^t×ev(is+1:ie,1:iwidth)を計算する。次に、ブロック形式の（１＋ＵＢＵ^t）のＢを求める。diag（w2）=-diag(a(is:is+blk-1,is:is+blk-1))、i2=blk、を計算する。ｗ２の部分に対応する係数αを格納する。なお、TRL（w2）はw2の下三角行列、diag(x)はｘの対角要素部分である。
【００８０】
ステップＳ１１８においては、i2<1か否かを判断する。ステップＳ１１８の判断がＮＯの場合には、ステップＳ１１９において、w2の上三角部分にハウスホルダベクトルの内積×αを格納し、i1=i2-1と設定する。ステップＳ１２０では、i1<1か否かを判断し、ＮＯの場合には、ステップＳ１２１において、w2(i1,i2)=w2(i1,i1)×（a(is+i2:n,is+i2-1）^t×a(is+i2:n,is+i1-1))、i1=i1-1を計算し、ステップＳ１２０に戻る。ステップＳ１２０の判断がＹＥＳの場合には、ステップＳ１２２において、i2=i2-1と設定し、ステップＳ１１８に戻る。ステップＳ１１８の判断がＹＥＳの場合には、ステップＳ１２３において、i1=blk-2とし、以下、２重ループで展開係数を計算する。三角行列の上辺を右から左に決定し、積み上げるように計算する。これは左からハウスホルダ変換をかけたときの展開を足して係数を決めることに相当する。ステップＳ１２４において、i1<1か否かを判断する。ステップＳ１２４の判断がＮＯの場合には、ステップＳ１２５においては、i2=blkと設定し、ステップＳ１２６においては、i2<i1+1か否かを判断する。ステップＳ１２６の判断がＮＯの場合には、ステップＳ１２７において、上辺の要素を左から右に決めていく。一つ前の係数を利用する。すなわち、w2(i1,i2)=w2(i1,i2)+w2(i1,i1+1:i2-1)×w2(i1+1:i2-1,i2)、i2=i2-1を計算し、ステップＳ１２６に戻る。ステップＳ１２６の判断がＹＥＳの場合には、ステップＳ１２８において、i1=i1-1と設定し、ステップＳ１２４に戻る。ステップＳ１２４の判断がＹＥＳの場合には、ステップＳ１２９に進み、i2=blkを設定し、ステップＳ１３０において、i2<1か否かを判断する。ステップＳ１３０の判断がＮＯの場合には、ステップＳ１３１において、以下のループで不足していた係数αをかける。まず、i1=i2-1と設定し、ステップＳ１３２において、i1<1か否かを判断する。ステップＳ１３２における判断がＮＯの場合には、ステップＳ１３３において、w2(i1,i2)=w2(i1,i2)×w2(i2,i2)、i1=i1-1を計算し、ステップＳ１３２に戻る。ステップＳ１３２の判断がＹＥＳの場合には、ステップＳ１３４において、i2=i2-1と設定し、ステップＳ１３０に戻る。ステップＳ１３０の判断がＹＥＳの場合、ステップＳ１３５において、ＢＵ^tを計算してＷに格納する。Ｗ（1:blk,1:iwidth）=TRU（w2）×Ｗ（1:blk,1:iwidth）を計算する。そして、Ｕの上部の三角行列と下部の長方形の部分とＷに格納したＢＵ^tを使って（１＋ＵＢＵ^t）×ＥＶを計算する。すなわち、ev(ie+1:n,1:width）=ev(ie+1:n,1:width)+a(ie+1:n,is:ie)×W(1:blk,1:width)、ev(is+1:ie,1:width)=ev(is+1:ie,1:width)+TRL（a(is+1:ie,is+1:ie))×W（1:blk-1,1:width)を計算し、ステップＳ１１５に戻る。
【００８１】
本発明によれば、共有メモリ型スカラ並列計算機を用いて高性能かつスケーラビリティのある固有値・固有ベクトルの並列演算方法を提供できる。
本発明の実施形態によれば、特に、固有ベクトルの変換に関する計算は従来の方法に比べて１０倍程度スカラマシンで高速。これらのアルゴリズムを使った実対称行列の固有値・固有ベクトルは、３重対角行列の固有値・固有ベクトルを同じスツルムの方法と逆反復法で求める方法を利用したもので、SUNの数値計算ライブラリSUN Performance libraryの機能に比べて、７個のＣＰＵで６．７倍高速。SUNの他のルーチン（機能が劣る：固有値・固有ベクトルを選択的には求められない）３重対角行列の固有値・固有ベクトルの計算をdivide ＆ conquer法で行うものと比べても２．３倍高速。
【００８２】
これらのアルゴリズムを使ったエルミート行列の固有値・固有ベクトルは、３重対角行列の固有値・固有ベクトルを同じスツルムの方法と逆反復法で求める方法を利用したもので、SUNの数値計算ライブラリSUN Performance libraryの機能に比べて７個のＣＰＵで４．８倍高速。SUNの他のルーチン（機能が劣る：固有値を選択的には求められない）３重対角行列の固有値・固有ベクトルの計算をdivide ＆ conquer法で行うものと比べても３．８倍高速。
【００８３】
なお、行列計算の基本的アルゴリズムに関しては、以下の教科書を参照されたい。
G. H. Golub and C. F. Van Loan "Matrix Computations" third edition The Johns Hopkins University Press, 1996
また、３重対角化の並列計算については、以下の文献を参照されたい。
J. Choi, J. J. Dongarra, and D. W. Walker, "THE DESIGN OF A PARALLEL DENSE LINEAR ALGEBRA SOFTWARE LIBRARY: REDUCTION TO HESSENBERG, TRIDIAGONAL, AND BIDIAGONAL FORM", Engineering Physics and Mathematics Division, Mathematical Sciences Section, prepared by the Oak Ridge National Laboratory managed by Martin Marietta Energy System, Inc., for the U.S. DEPARTMENT OF ENERGY under Contract No. DE-AC05-84OR21400, ORNL/TM-12472.（付記１）共有メモリ型スカラ並列計算機用固有値問題の並列処理方法を実現させるプログラムであって、
演算すべき実対称行列あるいはエルミート行列をブロック化し、ブロック化した部分をメモリの作業域にコピーして、ブロック間の積を利用して３重対角化するステップと、
３重対角化された行列から固有値と固有ベクトルを求めるステップと、
該３重対角化された行列から求められた固有ベクトルを、所定のブロック幅の行列演算の並列演算となるように変形されたハウスホルダー変換によって変換し、元の行列の固有ベクトルを求めるステップと、
を備えることを特徴とする並列処理方法を共有メモリ型スカラ並列計算機に実現させるプログラム。
【００８４】
（付記２）前記３重対角化ステップにおいて、ブロック化された部分の更新を再帰的プログラムで行うことを特徴とする付記１に記載のプログラム。
（付記３）前記３重対角化ステップにおいて、ブロック化された部分は、キャッシュメモリのページを跨らないように小さなブロックに分割して格納し、各プロセッサで並列に演算を行うことによって処理されることを特徴とする付記１に記載のプログラム。
【００８５】
（付記４）前記固有ベクトルを求めるステップにおいて、前記ハウスホルダー変換のための行列は、各プロセッサで処理できるような小さな係数行列となる上三角行列を各プロセッサで同時に生成されることを特徴とする付記１に記載のプログラム。
【００８６】
（付記５）前記固有ベクトルを求めるステップにおいて、前記固有ベクトルは、格納されている２次元配列の２次元目をプロセッサの数に対応して均等に分割し、おのおのの分割された領域を各プロセッサに割り当てることによって、演算されることを特徴とする付記１に記載のプログラム。
【００８７】
（付記６）共有メモリ型スカラ並列計算機用固有値問題の並列処理方法であって、
演算すべき実対称行列あるいはエルミート行列をブロック化し、ブロック化した部分をメモリの作業域にコピーして、ブロック間の積を利用して３重対角化するステップと、
３重対角化された行列から固有値と固有ベクトルを求めるステップと、
該３重対角化された行列から求められた固有ベクトルを、所定のブロック幅の行列演算の並列演算となるように変形されたハウスホルダー変換によって変換し、元の行列の固有ベクトルを求めるステップと、
を備えることを特徴とする並列処理方法。
【００８８】
（付記７）前記３重対角化ステップにおいて、ブロック化された部分の更新を再帰的プログラムで行うことを特徴とする付記６に記載の並列処理方法。
（付記８）前記３重対角化ステップにおいて、ブロック化された部分は、キャッシュメモリのページを跨らないように小さなブロックに分割して格納し、各プロセッサで並列に演算を行うことによって処理されることを特徴とする付記６に記載の並列処理方法。
【００８９】
（付記９）前記固有ベクトルを求めるステップにおいて、前記ハウスホルダー変換のための行列は、各プロセッサで処理できるような小さな係数行列となる上三角行列を各プロセッサで同時に生成されることを特徴とする付記６に記載の並列処理方法。
【００９０】
（付記１０）前記固有ベクトルを求めるステップにおいて、前記固有ベクトルは、格納されている２次元配列の２次元目をプロセッサの数に対応して均等に分割し、おのおのの分割された領域を各プロセッサに割り当てることによって、演算されることを特徴とする付記６に記載の並列処理方法。
【００９１】
【発明の効果】
高性能かつスケーラビリティのある固有値・固有ベクトルの解法を実現できる。
【図面の簡単な説明】
【図１】本発明の実施形態が前提とする共有メモリ型スカラ計算機のハードウェア構成を示した図である。
【図２】本発明の実施形態のアルゴリズムを説明する図（その１）である。
【図３】本発明の実施形態のアルゴリズムを説明する図（その２）である。
【図４】本発明の実施形態のアルゴリズムを説明する図（その３）である。
【図５】本発明の実施形態のアルゴリズムを説明する図（その４）である。
【図６】本発明の実施形態のアルゴリズムを説明する図（その５）である。
【図７】本発明の実施形態のアルゴリズムを説明する図（その６）である。
【図８】本発明の実施形態のアルゴリズムを説明する図（その７）である。
【図９】本発明の実施形態のアルゴリズムを説明する図（その８）である。
【図１０】本発明の実施形態のアルゴリズムを説明する図（その９）である。
【図１１】本発明の実施形態のアルゴリズムを説明する図（その１０）である。
【図１２】本発明の実施形態に従ったルーチンの疑似コード（その１）である。
【図１３】本発明の実施形態に従ったルーチンの疑似コード（その２）である。
【図１４】本発明の実施形態に従ったルーチンの疑似コード（その３）である。
【図１５】本発明の実施形態に従ったルーチンの疑似コード（その４）である。
【図１６】本発明の実施形態に従ったルーチンの疑似コード（その５）である。
【図１７】本発明の実施形態に従ったルーチンの疑似コード（その６）である。
【図１８】本発明の実施形態に従ったルーチンの疑似コード（その７）である。
【図１９】疑似コードの処理をフローチャートで示した図（その１）である。
【図２０】疑似コードの処理をフローチャートで示した図（その２）である。
【図２１】疑似コードの処理をフローチャートで示した図（その３）である。
【図２２】疑似コードの処理をフローチャートで示した図（その４）である。
【図２３】疑似コードの処理をフローチャートで示した図（その５）である。
【図２４】疑似コードの処理をフローチャートで示した図（その６）である。
【図２５】疑似コードの処理をフローチャートで示した図（その７）である。
【図２６】疑似コードの処理をフローチャートで示した図（その８）である。
【図２７】疑似コードの処理をフローチャートで示した図（その９）である。
【図２８】疑似コードの処理をフローチャートで示した図（その１０）である。
【図２９】疑似コードの処理をフローチャートで示した図（その１１）である。
【符号の説明】
１０−１〜１０−ｎ プロセッサ
１１−１〜１１−ｎ メモリモジュール
１２ 相互結合網
１３−１〜１３−ｎ ２次キャッシュメモリ

Claims (4)

1. 共有メモリ型スカラ並列計算機用固有値問題の並列処理方法を実現させるプログラムであって、
   演算すべき実対称行列あるいはエルミート行列Ａ（１：ｋ、１：ｋ）（ｋは、整数）の下三角部分を、列ブロックＡ（（ｉ−１）＊ｉｂｌｋ＋２：ｋ、（ｉ−１）＊ｉｂｌｋ＋１：ｉ＊ｉｂｌｋ）＝Ｕ、残りブロックＡ（（ｉ−１）＊ｉｂｌｋ＋２：ｋ、ｉ＊ｉｂｌｋ＋１：ｋ）＝Ａ _ｉ （ｉは、繰り返し数、ｉｂｌｋは、ブロック幅）とブロック化し、ブロック化した部分Ｕをメモリの作業域にコピーして、ブロック間の積を利用して３重対角化するステップと、
   ３重対角化された行列から固有値と固有ベクトルを求めるステップと、
   前記共有メモリ型並列計算機の各プロセッサで、該３重対角化された行列から求められた固有ベクトルＥＶ（１：ｎ、１＋（ｋ−１）＊ｎｂ：ｋ＊ｎｂ）（求める固有ベクトルをｎｅｖとしてｐ個のプロセッサで均等処理する個数をｎｂとする）に、所定のブロック幅ｎｂｌｋの行列演算をこれら分割された固有ベクトルに各々演算して並列計算するために、
   正方行列ＡのＡ（１：ｎ、１：ｎ−２）の下三角部分に実対称行列またはエルミート行列の３重対角化で用いたハウスホルダー変換１＋α _ｋ ＊ｕ _ｋ ｕ _ｋ ^ｔ （以下、エルミート行列の場合ｔ→Ｈ：Ｈは転置と複素共役）のα _ｋ を対角要素Ａ（ｋ、ｋ）に、また、ハウスホルダーベクトルｕ _ｋ をＡ（ｋ＋１：ｎ、ｋ）に格納しておき（１からｎ−２までを均等に分割したブロック幅がｎｂｌｋ）、１から始まる繰り返し数ｍに応じてＡ（１：ｎ、ｎ−２−ｍ＊ｎｂｌｋ＋１：ｎ−２−（ｍ−１）＊ｎｂｌｋ）部分の対角より下に格納されたハウスホルダーベクトルの束をＵ _ｍ とし、ｎｂｌｋ×ｎｂｌｋの行列Ｂを各プロセッサに取得し、Ａ（１：ｎ、ｎ−２−ｍ＊ｎｂｌｋ＋１：ｎ−２−（ｍ−１）＊ｎｂｌｋ）に格納されているハウスホルダーベクトルを小さい方から、１からｎｂｌｋまで番号を付け直した番号で示す添え字ｉ、ｊ（ｉ＜ｊ）で表す対応するハウスホルダーベクトルｕ _ｉ とｕ _ｊ の内積にα _ｉ を掛けたものを上三角部分Ｂ（ｉ、ｊ）に格納し、同じく付け直した番号ｊに対応するハウスホルダーベクトルに対応するα _ｊ をＢの対角要素Ｂ（ｉ、ｊ）に格納し、これらデータをさらに上三角部分で更新することで、各プロセッサで同じ計算を行い同じ行列Ｂを構成するステップと、
   各プロセッサで行列Ｂを利用して、（１＋Ｕ _ｍ ＢＵ _ｍ ^ｔ ）なる行列演算の形のハウスホルダー変換に変形し、固有ベクトルＥＶ（１：ｎ、１＋（ｋ−１）＊ｎｂ：ｋ＊ｎｂ）に左から行列演算の形で各プロセッサで並列演算を行い、これらの処理を繰り返し数ｍを＋１して繰り返すことで固有ベクトルを計算するステップと、
   を備えることを特徴とする並列処理方法を共有メモリ型スカラ並列計算機に実現させるプログラム。
2. 前記３重対角化ステップにおいて、ブロック化された部分の更新を再帰的プログラムで行うことを特徴とする請求項１に記載のプログラム。
3. 前記３重対角化ステップにおいて、ブロック化された部分は、キャッシュメモリのページを跨らないように小さなブロックに分割して格納し、各プロセッサで並列に演算を行うことによって処理されることを特徴とする請求項１に記載のプログラム。
4. 前記固有ベクトルを求めるステップにおいて、前記固有ベクトルは、格納されている２次元配列の２次元目をプロセッサの数に対応して均等に分割し、おのおのの分割された領域を各プロセッサに割り当てることによって、演算されることを特徴とする請求項１に記載のプログラム。

Images