JP4037303B2 - Parallel processing method of eigenvalue problem for shared memory type scalar parallel computer. - Google Patents

Description

ｂ）ハウスホルダー変換を構成する情報の格納
固有ベクトルを計算するときに３重対角化で使ったハウスホルダー変換が必要になる。このため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.（付記１）共有メモリ型スカラ並列計算機用固有値問題の並列処理方法を実現させるプログラムであって、
演算すべき実対称行列あるいはエルミート行列をブロック化し、ブロック化した部分をメモリの作業域にコピーして、ブロック間の積を利用して３重対角化するステップと、
３重対角化された行列から固有値と固有ベクトルを求めるステップと、
該３重対角化された行列から求められた固有ベクトルを、所定のブロック幅の行列演算の並列演算となるように変形されたハウスホルダー変換によって変換し、元の行列の固有ベクトルを求めるステップと、
を備えることを特徴とする並列処理方法を共有メモリ型スカラ並列計算機に実現させるプログラム。
【００８４】
（付記２）前記３重対角化ステップにおいて、ブロック化された部分の更新を再帰的プログラムで行うことを特徴とする付記１に記載のプログラム。
（付記３）前記３重対角化ステップにおいて、ブロック化された部分は、キャッシュメモリのページを跨らないように小さなブロックに分割して格納し、各プロセッサで並列に演算を行うことによって処理されることを特徴とする付記１に記載のプログラム。
【００８５】
（付記４）前記固有ベクトルを求めるステップにおいて、前記ハウスホルダー変換のための行列は、各プロセッサで処理できるような小さな係数行列となる上三角行列を各プロセッサで同時に生成されることを特徴とする付記１に記載のプログラム。
【００８６】
（付記５）前記固有ベクトルを求めるステップにおいて、前記固有ベクトルは、格納されている２次元配列の２次元目をプロセッサの数に対応して均等に分割し、おのおのの分割された領域を各プロセッサに割り当てることによって、演算されることを特徴とする付記１に記載のプログラム。
【００８７】
（付記６）共有メモリ型スカラ並列計算機用固有値問題の並列処理方法であって、
演算すべき実対称行列あるいはエルミート行列をブロック化し、ブロック化した部分をメモリの作業域にコピーして、ブロック間の積を利用して３重対角化するステップと、
３重対角化された行列から固有値と固有ベクトルを求めるステップと、
該３重対角化された行列から求められた固有ベクトルを、所定のブロック幅の行列演算の並列演算となるように変形されたハウスホルダー変換によって変換し、元の行列の固有ベクトルを求めるステップと、
を備えることを特徴とする並列処理方法。
【００８８】
（付記７）前記３重対角化ステップにおいて、ブロック化された部分の更新を再帰的プログラムで行うことを特徴とする付記６に記載の並列処理方法。
（付記８）前記３重対角化ステップにおいて、ブロック化された部分は、キャッシュメモリのページを跨らないように小さなブロックに分割して格納し、各プロセッサで並列に演算を行うことによって処理されることを特徴とする付記６に記載の並列処理方法。
【００８９】
（付記９）前記固有ベクトルを求めるステップにおいて、前記ハウスホルダー変換のための行列は、各プロセッサで処理できるような小さな係数行列となる上三角行列を各プロセッサで同時に生成されることを特徴とする付記６に記載の並列処理方法。
【００９０】
（付記１０）前記固有ベクトルを求めるステップにおいて、前記固有ベクトルは、格納されている２次元配列の２次元目をプロセッサの数に対応して均等に分割し、おのおのの分割された領域を各プロセッサに割り当てることによって、演算されることを特徴とする付記６に記載の並列処理方法。
【００９１】
【発明の効果】
高性能かつスケーラビリティのある固有値・固有ベクトルの解法を実現できる。
【図面の簡単な説明】
【図１】本発明の実施形態が前提とする共有メモリ型スカラ計算機のハードウェア構成を示した図である。
【図２】本発明の実施形態のアルゴリズムを説明する図（その１）である。
【図３】本発明の実施形態のアルゴリズムを説明する図（その２）である。
【図４】本発明の実施形態のアルゴリズムを説明する図（その３）である。
【図５】本発明の実施形態のアルゴリズムを説明する図（その４）である。
【図６】本発明の実施形態のアルゴリズムを説明する図（その５）である。
【図７】本発明の実施形態のアルゴリズムを説明する図（その６）である。
【図８】本発明の実施形態のアルゴリズムを説明する図（その７）である。
【図９】本発明の実施形態のアルゴリズムを説明する図（その８）である。
【図１０】本発明の実施形態のアルゴリズムを説明する図（その９）である。
【図１１】本発明の実施形態のアルゴリズムを説明する図（その１０）である。
【図１２】本発明の実施形態に従ったルーチンの疑似コード（その１）である。
【図１３】本発明の実施形態に従ったルーチンの疑似コード（その２）である。
【図１４】本発明の実施形態に従ったルーチンの疑似コード（その３）である。
【図１５】本発明の実施形態に従ったルーチンの疑似コード（その４）である。
【図１６】本発明の実施形態に従ったルーチンの疑似コード（その５）である。
【図１７】本発明の実施形態に従ったルーチンの疑似コード（その６）である。
【図１８】本発明の実施形態に従ったルーチンの疑似コード（その７）である。
【図１９】疑似コードの処理をフローチャートで示した図（その１）である。
【図２０】疑似コードの処理をフローチャートで示した図（その２）である。
【図２１】疑似コードの処理をフローチャートで示した図（その３）である。
【図２２】疑似コードの処理をフローチャートで示した図（その４）である。
【図２３】疑似コードの処理をフローチャートで示した図（その５）である。
【図２４】疑似コードの処理をフローチャートで示した図（その６）である。
【図２５】疑似コードの処理をフローチャートで示した図（その７）である。
【図２６】疑似コードの処理をフローチャートで示した図（その８）である。
【図２７】疑似コードの処理をフローチャートで示した図（その９）である。
【図２８】疑似コードの処理をフローチャートで示した図（その１０）である。
【図２９】疑似コードの処理をフローチャートで示した図（その１１）である。
【符号の説明】
１０−１〜１０−ｎ プロセッサ
１１−１〜１１−ｎ メモリモジュール
１２ 相互結合網
１３−１〜１３−ｎ ２次キャッシュメモリ[0001]
BACKGROUND OF THE INVENTION
The present invention relates to matrix operations in a shared memory scalar parallel computer.
[0002]
[Prior art]
Real symmetric matrices (matrix elements that are real numbers and remain unchanged even if transposed) and Hermitian matrices (matrix elements that are complex numbers and remain unchanged even if they are conjugate / transposed) for vector parallel computers The eigenvalue problem (when A is a matrix, λ is a constant, and I is a unit matrix, finding λ and its eigenvector where det | A−λI | = 0) is tridiagonalization (diagonal components and their First, do a matrix with only the components on both sides. The eigenvalue problem of this tridiagonal matrix is solved by the multisection method, eigenvalues are obtained, and then eigenvectors are obtained by the inverse iteration method. After that, householder transformation was performed on the eigenvectors to find the eigenvectors of the original eigenvalue problem.
[0003]
[Problems to be solved by the invention]
In a vector computer, the eigenvalue problem is calculated on the assumption that the memory access system operates fast, but in the case of a shared memory scalar computer, the number of times the shared memory is accessed as the calculation matrix increases. There is a problem that the performance of a computer is greatly impaired by accessing a shared memory with a slow access speed. Therefore, it is necessary to perform matrix calculation by effectively using a cache memory having a high access speed provided in each processor of the shared memory type scalar computer. In other words, if you calculate for each row or column of the matrix, access to the shared memory will increase, so after the matrix is blocked and the data stored in the cache memory in each processor is processed to the maximum, An algorithm localized to each processor is required that reduces the number of accesses to the shared memory by accessing the shared memory.
[0004]
That is, since a shared memory type parallel computer does not have a powerful memory access capability like a vector computer, it is necessary to devise an algorithm for increasing the amount of calculation for memory access.
[0005]
An object of the present invention is to provide a parallel processing method of operations that can perform an eigenvalue problem at high speed in a shared memory scalar computer.
[0006]
[Means for Solving the Problems]
The parallel processing method of the present invention is a program for realizing a parallel processing method for an eigenvalue problem for a shared memory type scalar parallel computer, which blocks a real symmetric matrix or Hermitian matrix to be operated, and the block portion is a memory operation. Copying to a region and tridiagonalizing using a product between blocks, obtaining eigenvalues and eigenvectors from the tridiagonalized matrix, and from the tridiagonalized matrix Transforming the obtained eigenvector by Householder transformation modified so as to be a parallel operation of a matrix operation having a predetermined block width to obtain an eigenvector of the original matrix.
[0007]
According to the present invention, since the eigenvalue problem can be solved by operations as localized as possible in each processor of the shared memory type scalar computer, the delay time caused by each processor frequently accessing the shared memory is minimized. It is possible to perform operations that effectively suppress the power of parallel computing.
[0008]
DETAILED DESCRIPTION OF THE INVENTION
In the embodiment of the present invention, the algorithm of the tridiagonalization part in the calculation of the eigenvalue problem employs a block algorithm, and recursively performs the algorithm of the part for obtaining the block part, and the update part Increase the computation density at. In addition, the matrix vector product portion is made to be a system that can use continuous access by utilizing symmetry so as to avoid memory access from accessing the pages of the cache memory step by step. When data is read across the pages of the cache memory, it is necessary to access the cache memory twice without being able to read the data at a time. Therefore, as described above, the page is not covered with data.
[0009]
In addition, the part that performs householder transformation on the eigenvectors of the tridiagonalized matrix and obtains the eigenvectors of the original matrix bundles the householder transformations about every 80 times and performs an operation consisting of three matrix elements. This increases the operation density.
[0010]
In the embodiment of the present invention, the method for obtaining eigenvalues from a tridiagonalized matrix and the method for obtaining eigenvectors of a tridiagonalized matrix are based on conventional methods.
[0011]
FIG. 1 is a diagram showing a hardware configuration of a shared memory scalar computer on which the embodiment of the present invention is premised.
The processors 10-1 to 10-n have a primary cache memory, and this primary cache memory may be incorporated in the processor. The processors 10-1 to 10-n are provided with secondary cache memories 13-1 to 13-n, and the secondary cache memories 13-1 to 13-n are coupled to the interconnection network 12. . Further, the mutual coupling network 12 is provided with memory modules 11-1 to 11-n which are shared memories, and the processors 10-1 to 10-n read data necessary for calculation from here, and the mutual coupling network 12 The secondary cache memories 13-1 to 13-n or the primary cache memory are used for calculation.
[0012]
In this case, data is read from the memory modules 11-1 to 11-n to the secondary cache memories 13-1 to 13-n or the primary cache memory, the secondary cache memories 13-1 to 13-n, or 1 Writing the calculated data from the next cache memory to the memory modules 11-1 to 11-n is much slower than the calculation speed of the processors 10-1 to 10-n. Therefore, if such writing and reading frequently occur, the performance of the entire computer is degraded.
[0013]
Therefore, in order to keep the performance of the entire computer high, access to the memory modules 11-1 to 11-n is reduced as much as possible, the secondary cache memories 13-1 to 13-n, the primary cache memory, and the processor 10 An algorithm that performs as many operations as possible in a local system consisting of -1 to 10-n is required.
<Method for obtaining eigenvalues and eigenvectors>
1.3 Double diagonalization
1) Tridiagonalization
a) Mathematical algorithm for blocked tridiagonalization
The matrix is tridiagonalized for each block width. The matrix is divided into blocks and tridiagonalization is performed on a block basis by the following algorithm.
[0014]
2-11 is a figure explaining the algorithm of embodiment of this invention.
FIG. 2 is a diagram when the m-th block is processed. The block is a column and a rectangle in the row direction indicated by dotted lines in FIG.
[0015]
For the last block, all processing is completed by applying the algorithm to the left-justified block with the block width −2 as the block width and leaving the 2 × 2 matrix in the lower right.
[0016]
do i = 1, blks
step1: A_nMake a house holder vector u from n + 1 row vectors.
step2: v_i= A_{n + i}u, w_i= V_i-u (u^tv) Calculate / 2.
step3: U_i= (U_i-1, U_i), W_i= (W_i-1, W_i) And update. (Where (U_i-1, U_i) Is the matrix U_i-1Add one column and matrix U_iThus, the matrix is expanded by one column. )
step4: if (i <blks) then
A_nUpdate the n + i + 1 column.
[0017]

b) Storing the information that makes up the householder transformation
When calculating eigenvectors, the householder transformation used in tridiagonalization is required. For this reason U_nAnd α are stored in the positions of vectors constituting the householder transformation. α is stored at the position of the corresponding diagonal element.
c) U_iTo find out efficiently
It is necessary to update the next vector for householder transformation in order to tridiagonalize the block part. In order to perform these calculations as locally as possible, the block width portion is copied to the work area, tridiagonalized, and stored in the original area. Instead of updating the next column vector each time, the calculation is performed in the form of a matrix product and the operation density is increased. For this reason, tridiagonalization of the block portion is performed by recursive programming.
[0018]
recursive subroutine trid (width, block area pointer)
if (width <10) then
c Tridiagonalization of the block portion of the width.
[0019]
From vector u and matrix vector product required for householder transformation, v_iAnd w_iCreate
U, W to u_i, W_iJoin.
[0020]
else
c Divide the block width in half
c Tridiagonalization is performed on the first half block.
[0021]
call trid (first half width, first half area)
c Divide the block and update the second half from the dividing line
B = B-UW^t-WU^tUpdate with.
c Then, the latter half is tridiagonalized.
[0022]
call trid (second half width, second half area)
return
end
As shown in FIG. 3, the block portion is copied to the work area U, and the block portion is tridiagonalized by a recursive program. Because of the recursive program, the first half of FIG. 3 is tridiagonalized when the recursive program is called for the first half update process. The second half is updated by the first half and then tridiagonalized.
[0023]
As shown in FIG. 4, when the recursive program is called up to a depth of 2, the hatched portion in FIG. 4A is updated as B in the first first half process, and then the shaded portion in FIG. Finally, FIG. 4E is updated. In the parallelization at the time of update, the block matrix of the update part is equally divided in the row vector direction, and the update of each part is performed in parallel by a plurality of processors.
After the calculation of FIG. 4 (a), the calculation of FIG. 4 (b) is performed, after the calculation of FIG. 4 (c), the calculation of FIG. 4 (d) is performed, and after the calculation of FIG. 4 (f) is calculated.
[0024]
As shown in FIG. 5, when updating the hatched portion of U, the horizontal line of u and the vertical line portion of W are referred to. This can increase the calculation density. That is,
A_{n + k}= A_n-U_kW_k ^t-W_kU_k ^t   ....... (**)
V_nCan be calculated as follows: At this time, the reference pattern of U and W is
v_k= A_nu_k-U_k-1W_k-1 ^tu_k-W_k-1U_k-1 ^tu_k    ..... (***)
After updating U in FIGS. 5 (a) and (b), (c) and (d), (e) and (f), v_kIs calculated by a matrix vector product with reference to U and W. It is only a reference, and there is a part in which the update and reference of U are in common, and the reference can be made efficiently. A_nIn this way, the necessary parts are updated using U and W, and the performance density is improved by using (**). (***) is an extra operation, but can be ignored as long as the block width is reduced.
[0025]
For example, when performing parallel processing with four CPUs, (***) W_k-1 ^tu_k, U_k-1 ^tu_kIn this calculation, parallel calculation is performed by dividing the vertical line, and a product with the result is divided in the direction of the broken line and calculated in parallel.
v_i= A_nu_iParallel computation
As shown in FIG._nSymmetry, that is, A_n= A_n ^tDivide the second dimension with each processor using v_iPart A_n(*, Ns: ne)^tu_iIs calculated for each processor.
2) Parallelization on a shared memory scalar parallel computer
a) The storage areas for U and W are secured in the shared memory. Further, the block area to be tridiagonalized is copied to a separately secured work area, and tridiagonalization is performed in that area.
[0026]
In the above recursive program, the parts to be parallelized are the following parts:
(1) u needed to build householder transformation_iStep4 to ask for
A_n(*, n + i + 1) = A_n(*, n + i + 1) -U_iW_i(n + i + 1, *)^t-W_iU_i(n + i + 1, *)^t
Use to find the required vector
(2) v at step2_iThe part to calculate. This part is added to the following expression of (**) u_iTo calculate.
A_{n + k}= A_n-U_kW_k ^t-W_kU_k ^t
In this calculation, A_nAnd u_iProduct with and U_kW_k ^t-W_kU_k ^tAnd u_iParallelize the product with.
[0027]
Block part is copied to work area A_nThe necessary part of is not updated. A_nUsing the symmetry of, divide into long matrices in the column vector direction and calculate in parallel.
(3) In recursive program, block area update is
A_{n + k}= A_n-U_kW_k ^t-W_kU_k ^t
Update using. This reduces the calculation of (1).
3) Update part of step5
The update uses symmetry to calculate the lower half of the diagonal component. In the parallelization, when the number of CPUs is #CPU so that the load is balanced, the second dimension of the subarray in which the submatrix to be updated is stored is equally divided into 2 × # CPUs, and 1 to 2 × # Number with CPU. For the CPUs 1 to #CPU, the i-th part and the 2 × # CPU + 1−i-th part are updated in parallel by the i-th processor.
[0028]
Then copy to the upper half. This is similarly divided to balance the load. At this time, the portions other than the diagonal block portion are divided into small blocks to be copied so as not to cross the cache memory page.
A_{n + k}= A_n-U_kW_k ^t-W_kU_k ^tThe part to be updated is a lower triangular matrix part, which is divided into column blocks divided by the number of CPUs × 2 as shown in FIG. 7, and a pair is formed from both ends, which is updated by each CPU. FIG. 7 describes a case where four CPUs are provided.
[0029]
When the lower triangular part can be updated, 1-8 are transposed in the same pair to the upper triangular part and copied from u1 to u8.
At this time, the inside is divided into small square blocks and transposed using a cache. Similar to the update, the CPU performs parallel processing.
Explanation of the performance improvement by transposition on the cache
As shown in FIG. 8, transpose conversion is performed from 1 to the numerical order. For each square area, in the case of 1, the lower triangular part is copied to a continuous area of the memory, and it is accessed (transposed) in the row direction and stored in the upper triangular part of 1. 2 to 8 are copied into a square and stored in the storage area in the column direction (vertical direction) with reference to the row direction.
Calculation of two eigenvectors
a) Basic algorithm
Vector u_nAnd later (1-2 * uu^t/ (u^tu)) and multiply by vector.
[0030]
When tridiagonalization is performed, the original eigenvalue problem is as follows.
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 from the eigenvector y obtained by solving the tridiagonalized eigenvalue problem₁ ^tQ₂ ^t..... Q_n-3 ^tQ_n-2 ^tIt can be converted by calculating y.
b) Parallel computation of block algorithm and eigenvector transformation of embodiment of the present invention
When obtaining a large number or all of the eigenvectors, the eigenvectors of the tridiagonal matrix are equally associated with each CPU, and the above-described conversion is performed in parallel on each CPU. At this time, about 80 conversion matrices are collectively converted.
[0031]
Each transformation matrix Q_i ^tIs 1 + α_iu_iu_i ^tCan be written. When these products are expanded, the following form is obtained.
b_{i, j}Is the leftmost and rightmost u_iu_jThis is a summary of scalar coefficients excluding.
[0032]
[Expression 1]

[0033]
It becomes the form. b_{i, j}Becomes an upper triangular matrix. 1 + UBU for this^tIt becomes the form. By transforming into such a shape, the calculation density is improved to increase the speed. FIG. 9 is a schematic diagram of the matrix B.
[0034]
Although the above method has three steps, there are two matrices, U and B, because of the memory access, and B can be made a relatively small area, which is efficient.
up to m-1st b_{i, j}(1 + α_mU_mU_m ^t)
[0035]
[Expression 2]

[0036]
Replace the sum of the last terms with i and j
[0037]
[Equation 3]

[0038]
And can be changed.
B in the innermost bracket_{m, j}(J = m + 1,..., N + k). b_{m, m}Is α_mIt is.
Prepare a square work array W2 with a block width.₁U_iU_j ^tIs stored in the upper triangular part of w2 (i, j). Α for diagonal elements_iIs stored.
[0039]
Next, the above argument can be calculated in the manner of adding one line at a time in order from the lower right 2 × 2 upper triangular matrix.
If calculation is performed one by one from the rightmost element of the row, B is an upper triangular matrix, so there is no reference to the updated part, so calculation can be performed in the same region. In this way, a coefficient matrix located in the middle of the three matrix products can be calculated using only a very small area.
[0040]
FIG. 10 is a diagram schematically showing a method for obtaining the eigenvalue.
Let the block width be nbs.
First, the inner product α_ju_i・ U_jIs calculated and stored.
[0041]
Α for diagonal elements_iIs stored.
Next, the calculation is performed as follows.
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
FIG. 11 is a diagram schematically showing how the eigenvectors obtained above are converted into eigenvectors of the original matrix.
[0042]
The eigenvector is converted by the householder vector stored in the array A. The result is stored in W, with the hatched portion in FIG. 11 applied to the hatched portion of EV. Also, W2 is constructed from the A block matrix. The hatched portion of EV is updated from the product of W2 and W multiplied by the block portion of A.
Eigenvalues and eigenvectors of 3 Hermitian matrices
The algorithm for obtaining eigenvalues and eigenvectors of the Hermitian matrix replaces the transposition in the tridiagonalization of the real symmetric matrix with transposition + complex conjugate (t → H). Also, how to make a householder vector is to add the vector size so that it is a scalar multiple of the first element.
[0043]
The resulting tridiagonal matrix is a Hermitian matrix, which is scaled by a diagonal matrix of absolute value 1.
The diagonal matrix is created as follows.
d_i= 1.0, d_{i + 1}= H_{i + 1, i}/ | H_{i + 1}| * D_i
12-18 are pseudo code for routines according to an embodiment of the present invention.
[0044]
FIG. 12 is a subroutine for tridiagonalizing a real symmetric matrix.
The lower triangular part of the real symmetric matrix is stored in the array a. A tridiagonal matrix and a sub-diagonal part are stored in daig and sdiag. Information necessary for conversion is stored in the lower triangular part of a as an output.
[0045]
U stores a block to be tripled. V is an area for storing W.
nb is the number of blocks, and nbase indicates the start position of the block.
After executing the subroutine copy, the block to be tridiagonalized is stored in u (nbase + 1: n, 1: iblk) and the blktrid routine is called to perform LU decomposition. Then, u (nbase + 1: n, 1: iblk) that has been processed is written back to the original matrix a. In the lower processing, the last remaining block is tridiagonalized using a subroutine blktrid.
[0046]
FIG. 13 is a pseudo code of a subroutine that performs tridiagonalization.
This subroutine is a routine for tridiagonalizing the block matrix and is called recursively. nbase is an offset indicating the position of the block, and istart is an offset within the block of the reduced subblock to be recursively called, and indicates the position of the target subblock when it is recursively called at first. nwidth is the size of the sub-block.
[0047]
If nwidth is less than 10, call subroutine bunit; otherwise, store istart in istart2, store half of nwidth in nwidth2, tridiagonalize by subroutine blktrid, Take barrier synchronization.
[0048]
Further, the sum of istart and nwidth / 2 is stored in istart3, nwidth-nwidth / 2 is stored in nwidth3, and a value is given to is2, is3, ie2, ie3, is, and ie indicating the beginning and end of the block, respectively. Set, set len, iptr, calculate according to the formula in the figure, store in u (is: ie, is3: ie3), take barrier synchronization, call tridiagonalizing subroutine blktrid and process To finish the subroutine processing.
[0049]
FIG. 14 is a pseudo code of the internal routine of the tridiagonalization subroutine.
In the tridiagonalization internal subroutine bunit, after storing necessary information, the start iptr2 of the block, the width len, the start position is, and the end position ie are determined, and barrier synchronization is established. And u (is: ie, i)^t* Store (is: ie, i) in tmp and take barrier synchronization. Then, each value is calculated and stored in each array. Here, sum means sum and sqrt means square root. Finally, barrier synchronization is taken.
[0050]
Next, v (is: ie, i) is obtained, barrier synchronization is performed, len2, isx, iex, u, v are updated, and barrier synchronization is performed. Furthermore, v (is: ie, i) is updated to obtain barrier synchronization. Furthermore, v (is: ie, i)^t* Calculate u (is: ie, i), store it in tmp, and take barrier synchronization.
[0051]
Then, a value is set in beta, barrier synchronization is performed, v is updated by an operation using beta, and barrier synchronization is performed.
Then, when i <iblk and ptr2 <n−2, u (is: ie, i + 1) is updated. In other cases, u (is: ie, i; 1: i + 2) is changed by another formula. Update and end the process. Note that after this subroutine is executed, the thread assignment is released.
[0052]
FIG. 15 shows pseudo codes of a routine for updating the lower half of the matrix from u and v, a routine for updating the diagonal matrix portion, and a copy routine.
Here, nbase is an offset indicating a block position, and nwidth is a block width.
[0053]
In this subroutine update, after ensuring the arrays a, u, v, barrier synchronization is established. Then, after setting blk, nbase2, len, is1, ee1, nbase3, isr, ier, a (ie1: n, is1: ie1) and a (ier + 1: n, isr: ier) are updated respectively. The subroutine trapdate is called twice, the barrier synchronization is taken, and the process is returned.
[0054]
In the subroutine trapdate, the block width for diagonal block update is stored in blk2, and an array of a, u, and v is secured.
is2 and ie2 are set, A (is2: ie, is2: ie2) is updated, and the process is returned.
[0055]
In the subroutine copy, len, is1, len1, nbase, isr, and lenr are set, bandcp is executed twice, and the process is returned.
FIG. 16 is a pseudo code of a routine for copying the updated lower triangular portion to the upper triangular portion.
[0056]
In the subroutine bandcp, nb, w, nn, and loopx are set, and TRL (a (is2: is2 + nnx-1, is2: is2 + nnx)) is set to TRL (w (1: nnx, 1: nnx)) in the do loop. And TRL (a (is2: is2 + nnx-1, is: is + nnx)) to TRL (w (i: nnx, 1: nnx))^tIs stored. Here, TRL indicates the lower triangular portion, and TRU indicates the upper triangular portion.
Then, w (1: nnx, 1: nnx) and a (is2: is2 + nnx, is3: is3 + nnx-1) are updated, and then w (1: ny, 1: nx) and a (is2: is2 + nnx, is3 : Update n).
[0057]
When the do loop is exited, the process is returned.
FIG. 17 is a pseudo code of a routine for converting the eigenvector of the tridiagonal matrix into the eigenvector of the original matrix.
[0058]
Here, the eigenvectors of the tridiagonal matrix are stored in ev (1: n, 1: nev). Further, a is an output of tridiagonalization, and information necessary for conversion is stored in the lower triangular part.
[0059]
Subroutine convev takes a and ev as an array.
Configure threads and perform parallel processing.
Barrier synchronization is established between threads, len, is, ie, and nevthrd are set, the convthrrd routine is called, and after returning, the barrier synchronization is established and the process is terminated.
[0060]
FIG. 18 is a pseudo code of a routine for converting eigenvectors.
In the subroutine convexthrd, the block width is stored in blk, and a, ev, w, and w2 are taken as an array.
[0061]
First, if the width is smaller than 0, the process returns to the original routine without processing. numblk and nfbs are set, and alpha (-a (i, i)) obtained by reversing the sign of the value stored in the diagonal element at the time of tridiagonalization is entered. Also, ev (i + 1: n, 1: iwidth) in x (1: iwidth)^t* A (i + 1: n, i) is input, and ev is updated using this and alpha and a. Furthermore, in the next do statement, is and ie are set, and w (1: blk, iwidth) is a (is + 1: n, is: ie)^t* Overwrite with ev (is + 1: n, 1: iwidth), and also change w (1: blk, 1: iwidth) to TRL (a (ie + 1: is, is: ie))^t* Update with ev (ie + 1: is, 1: iwidth). Here, TRL is the lower triangular matrix part.
[0062]
Also, the diagonal element vector of a (is: ie, is: ie) is stored in the diagonal element vector DIAG (w2) of w2.
In the next do statement, w2 (i1, i2) is changed to w2 (i1, i2) * (a (is + i2: n, is + i2-1)^t* A (is + i2: n, is + i1-1)). In the next do statement, w2 (i1, i2) is updated with w2 (i1, i2) + w2 (i1, i1 + 1: i2-1) * w2 (i1 + 1: i2-1, i2).
[0063]
Furthermore, in the next do statement, w2 (i1, i2) is updated with w2 (i1, i2) * w2 (i2, i2), and w (1: blk, 1: iwidth), ev (is + n: n, 1 : Iwidth), ev (ie + 1: is, 1: iwidth), and returns the process to the original routine.
[0064]
FIGS. 19 to 29 are flowcharts showing the pseudo code processing.
FIG. 19 is a flowchart of a subroutine trid for tridiagonalizing a real symmetric matrix. In step S10, the shared array A (k, n), diag (n) and sdiag (n) are input as a subroutine. diag and sdiag return the diagonal elements and sub-diagonal elements of the calculated tridiagonal matrix as output. The work areas U (n + 1, iblk) and v (n + 1, iblk) are secured inside the routine and used with the shared attribute. In step S11, a thread is generated, and the thread number assigned to each thread is set to the total thread number nothrd in the local area numthr in each thread, and the following is set in each thread. Set the block width to iblk and set nb = (n−2 + iblk−1) / iblk, nbase = 0, i = 1. In step S12, it is determined whether i> nb-1. If the determination in step S12 is yes, the process proceeds to step S19. If the determination in step S12 is NO, in step S13, nbase = (i−1) × iblk, istart = 1, and nwidth = iblk are set. In step S14, the subroutine copy is called to copy the lower triangular portion to the upper triangular portion. In step S15, the block tridiagonalization target area is copied to the work area U. That is, U (nbase + 1: n, 1: iblk) ← A (nbase + 1: n, nbase + 1: nbase + iblk) is executed. In step S16, the subroutine blktrid is called to tridiagonalize the portion copied to U. (Istart = 1, block width passes iblk). In step S17, the tridiagonalized one is returned to the sequence A. That is, A (nbase + 1: n, nbase + 1: nbase + iblk) ← U (nbase + 1: n, 1: iblk) is performed. In step S18, the subroutine update is called to update the lower triangular part of A (nbase + iblk: n, nbase + iblk: n), and the process returns to step S12.
[0065]
In step S19, nbase = (nb-1) × iblk, istart = 1, and iblk2 = n-nbase are set. In step S20, the block tridiagonalization target area is copied to the work area U. That is, U (nbase + 1: n, 1: nwidth) ← A (nbase + 1: n, nbase + 1: n) is performed. In step S21, the subroutine blktrid is called to tridiagonalize the portion copied to U (istart = 1, the block width is iblk2). In step S22, the tridiagonalized one is returned to the array A. That is, A (nbase + 1: n, nbase + 1: n) ← U (nbase + 1: n, 1: nwidth) is performed. In step S23, the thread generated for parallel processing is deleted, and the subroutine is terminated.
[0066]
FIG. 20 is a flowchart of subroutine blktrid. This subroutine is a recursive program.
This subroutine is called by the following statement:
subroutine blktrid (A, k, n, dig, sdig, nbase, istart, nwidth, U, V, nothrd, numthrd)
Here, nbase is an offset indicating the position of the block. istart is an offset within the block of the reduced block that is symmetric in the recursive call, and is initially 1 to indicate the position of the target block when it is recursively called. nwidth represents the block width. In step S25, it is determined whether nwidth <10. If the determination in step S25 is no, the process proceeds to step S27. If the determination in step S25 is yes, in step S26, the subroutine btunit is called to perform tridiagonalization, and the subroutine is exited. In step S27, iupdate2 = istart and nwidth2 = nwidth / 2 to be called are set and transferred while changing the target update position and block width for the recursive call. Pass the start position and block width of the reduced block. In step S28, the subroutine blktrid is recursively called. In step S29, barrier synchronization is established between threads. In step S30, the start point (is2, is3) and the end point (ie2, ie3) shared by each thread in the update are calculated. Istart3 = istart + nwidth / 2, nwidth3 = nwidth-nwidth / 2, is2 = istart2, ie2 = istart + nwidth2-1, is3 = istart3, ie3 = istart3 + nwidth3-1, iptr = nbase + istart3, len = Calculate (n-iptr + numthrd-1) / numthrd, is = iptr + (nothrd-1) * len + 1, ie = min (n, iptr + nothrd * len). In step S31, U (is: ie, is3: ie3) = U (is: ie, is3: ie3) −U (is: ie, is2: ie2) × W (is3: ie3, is2: ie2)^t-W (is: ie, is2: ie2) x U (is3: ie3, is2: ie2)^tIn step S32, barrier synchronization is established between threads. In step S33, the subroutine blktrid is recursively called to exit the subroutine.
[0067]
21 and 22 are flowcharts of a subroutine btunit which is an internal routine of the subroutine blktrid.
In step S35, tmp (numthrd), sigma, and alpha are assigned with the shared attribute. In step S36, it is determined whether nbase + istart> n-2. If the determination in step S36 is yes, the subroutine is exited. If the determination in step S36 is no, the process proceeds to step S38. In step S38, i = istart is set, and in step S39, it is determined whether i <= istart-1 + nwidth. If the determination in step S39 is no, the subroutine is exited. If the determination in step S39 is yes, in step S40, the start point is and end point ie shared by each thread are calculated. iptr2 = nbase + i, len = (n-iptr2 + numthrd-1) / numthrd, is = iptr2 + (nothrd-1) × len + 1, ie = min (n, iptr2 + nothrd × len). In step S41, barrier synchronization is taken, and in step S42, tmp (nothrd) = U (is: ie, i)^tXU (is: ie, i) is calculated, and barrier synchronization is established in step S43. In step S44, it is determined whether nothrd = 1. If the determination in step S44 is no, the process proceeds to step S46. If the determination in step S44 is yes, in step S45, the square root of the sum partially calculated in each thread is taken and calculation for tridiagonalization is performed (creation of a house holder vector). sigma = sqrt (sum (tmp (1: numthrd))), where SUM is the sum and sqrt is the square root. diag (iptr2) = u (iptr2, i), sdiag (iptr2) =-sigma, U (nbase + i + 1, i) = U (nbase + i + 1, i) + sign (u (nbase + i + 1, i) × sigma, alpha = 1.0 / (sigma × u (nbase + i + 1, i), U (iptr2, i) = alpha, and the process proceeds to step S46, where barrier synchronization is performed. In step S47, iptr3 = iptr2 + 1, and in step S48, V (is: ie, i) = A (iptr3: n, iptr2 + is: iptr2 + ie)^tU (ptr3: n, i) is calculated, and barrier synchronization is established in step S49.
[0068]
In step S50, V (is: ie, i) = alpha * V (is: ie, i) -V (is: ie, 1: i-1) * (U (iptr3: n, 1: i-1)^t× U (iptr3: n, i))-U (is: ie, 1: i-1) × (V (iptr3: n, 1: i-1)^t× U (iptr3: n, i)) is calculated, and barrier synchronization is established in step S51. In step S52, tmp (nothrd) = V (is: ie, i)^tX U (is: ie, i) is calculated, and barrier synchronization is established in step S53. In step S54, it is determined whether nothrd = 1. If the determination in step S54 is NO, the process proceeds to step S56. If the determination in step S54 is YES, the process proceeds to step S55. In step S55, beta = 0.5 * alpha * sum (tmp (1: numthrd)) is calculated. Here, sum is a symbol that takes the sum of vectors. In step S56, barrier synchronization is taken. In step S57, V (is: ie, i) = V (is: ie, i) −beta × U (is: ie, i) is calculated. In step S58, barrier synchronization is established. In step S59, it is determined whether or not ptr2 <n-2. If the determination is YES, U (is: ie, i + 1) = U (is: ie, i + 1) in step S60. -U (is: ie, istart: i) × V (i + 1, istart: 1)^t-V (is: ie, istart: i) × U (n + 1, istart: 1)^tAnd returns to step S39. If the determination in step S59 is NO, in step S61, U (is: ie, i + 1: i + 2) = U (is: ie, i + 1: I + 2) −U (is: ie , istart: i) x V (i + 1: n, istart: i)^t-V (is: ie, istart: i) x U (n + 1: n, istart: i)^tAnd exit the subroutine.
[0069]
FIG. 23 is a flowchart of the subroutine update.
In step S65, barrier synchronization is taken, and in step S66, a pair is formed by each thread, and the start point and end point for sharing the update are determined. That is, 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, and ier = min (n, isr + len−1). In step S67, A (ie1 + 1: n, is1: ie1) = A (ie1 + 1: n, is1 + 1: n, is1: ie1) −W (ie1 + 1: n, 1: blk) × U (Is1: ie1,1: blk)^t-U (ie1 + 1: n, 1: blk) x W (is1: ie1,1: blk)^t, A (ier + 1: n, isr: ier) = A (ier + 1: n, isr: ier) -W (ier + 1: n, 1: blk) × U (isr: ier, 1: blk)^t-U (ier + 1: n, 1: blk) x W (isr: ier, 1: blk)^tIn step S68, the subroutine trupdate is called to update the left half diagonal matrix portion. Pass is1, ie1, A, W, U. In step S69, the subroutine trupdate is called to update the right half diagonal matrix portion. Pass isr, ier, A, W, U. In step S70, barrier synchronization is taken and the subroutine is exited.
[0070]
FIG. 24 is a flowchart of subroutine trupdate (update of the diagonal matrix portion).
The update start point is and the end point ie are input. This is the start and end point used to update the rectangular part below the diagonal block before the subroutine is called.
[0071]
In step S75, the diagonal block update block width is set in blk2, and i = is is set. In step S76, it is determined whether i> ie-1. If the determination in step S76 is yes, the subroutine is exited. If the determination in step S76 is no, in step S77, the update start and end points for each thread are determined. That is, is2 = i, ie2 = min (i + blk2-1, ie-1), A (is2: ie-1, is2, ie2) = A (is2: ie-1, is2, ie2) -U (is2 : ie-1,1: blk) x W (is2: ie2,1: blk)^t-W (is2: ie-1,1: blk) x U (is2: ie2,1: blk)^tIn step S78, i = i + blk2 is set, and the process returns to step S76.
[0072]
FIG. 25 is a flowchart of subroutine copy.
In step S80, the starting point and width when copying is executed in parallel by making a pair in each thread are calculated. That is, 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)). In step S81, the subroutine bandcp is called. Copies the area created by the starting point is1 and width len1 on the left side of the pair. In step S82, a subroutine bandcp is called. Copies the area formed by the starting point isr and width lenr on the right side of the pair.
[0073]
FIG. 26 is a flowchart of the subroutine bandcp.
This routine performs copying while performing transposition on the cache using the small work area WX. The work area is WX (nb, nb), the start point is is, and the width is len.
[0074]
In step S85, nn = min (nb, len), loopx = (len + nn-1) / nn, and j = 1 are calculated. In step S86, it is determined whether j> loopx. If the determination in step S86 is yes, the subroutine ends. If the determination in step S86 is no, in step S87, the size of the diagonal block nnx to be copied to WX and its offset ip are determined. 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 (A (ip: ip + nnx-1, ip: ip + nnx-1)), TRU (A (ip: ip + nnx- 1, ip: ip + nnx-1)) = TRL (WX (1: nnx, i: nnx)), i = 1, is2 = ip, is3 = ip + nnx. However, TRU is the upper triangle and TRL is the lower triangle.
[0075]
In step S88, it is determined whether i> loopy-1. If the determination in step S88 is no, in step S89, the portion of nn × nnx is transposed and copied. That is, WX (1: nn, 1: nnx) = A (is3: is3 + nn-1, is2: is2 + nnx-1), A (is2: is2 + nnx-1, is3: is3 + nn-1) = WX (1, nn: 1, nnx)^t, Is3 = is3 + nn, and the process returns to step S88. If the determination in step S88 is yes, the last part is copied in step S90. That is, nn = n-is3 + 1, WX (1: nn, 1: nx) = A (is3: n, is2: is2 + nnx-1), A (is2: is2 + nnx-1, is3: n) = WX (1: nn, 1: nx) is calculated, and the process returns to step S86.
[0076]
FIG. 27 is a flowchart of the subroutine convev.
In this routine, the number nev of eigenvectors to be obtained and the house holder vector are stored in the lower half of a. The eigenvector of the tridiagonal matrix is stored in ev (k, nev).
[0077]
In step S95, a thread is generated. The total number of threads and the thread number (1 to numthrd) are set in numthr and nothrd in the local area of each thread. In step S96, barrier synchronization is taken. In step S97, a start point and an end point for sharing calculation are determined for each thread. That is, len = (nev + numthrd−1) / numthrd, is = (nothrd−1) × len + 1, ie = min (nev, nothrd × len), and width = ie-is + 1 are calculated. In step S98, the subroutine convexthrd is called to convert the tridiagonal matrix into the eigenvector of the original matrix. The area where the eigenvectors in charge of each thread are stored and the number of eigenvectors width are passed. In step S99, barrier synchronization is performed. In step S100, the generated thread is deleted, and the subroutine is exited.
[0078]
28 and 29 are flowcharts of the subroutine convexthrd.
This routine transforms the eigenvectors of the tridiagonal matrix assigned to each thread into the eigenvectors of the original matrix. Vectors and coefficients for restoring the householder transformation are stored in array A.
[0079]
In step S110, the block width is set to blk. The block width is about 80. In step S111, it is determined whether iwidth <0. If the determination in step S111 is yes, the subroutine is exited. If the determination in step S111 is no, the process proceeds to step S112. In step S112, the first block of transformation in the following loop is (1 + αuu).^t) Are calculated sequentially. First, numblk = (n−2 + blk−1) / blk and nfbs = n-2-blk × (numblk−1) are obtained. In step S113, it is determined whether i <n-2-nfbs + 1. If the determination in step S113 is YES, the process proceeds to step S114, where 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)^tAnd returns to step S113. If the determination in step S113 is NO, i = 1 is set in step S115, and it is determined in step S116 whether i> numblk-1. If the determination in step S116 is no, the subroutine is exited. If the determination in step S116 is YES, in block S117, the block format (1 + UBU^t) U^t× EV to U^tThe calculation is divided into the upper triangular matrix part at the left end and the rectangular part on the right side. That is, 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))^tXev (is + 1: ie, 1: iwidth) is calculated. Next, the block format (1 + UBU^tB) is obtained. diag (w2) =-diag (a (is: is + blk-1, is: is + blk-1)), i2 = blk is calculated. The coefficient α corresponding to the portion w2 is stored. TRL (w2) is a lower triangular matrix of w2, and diag (x) is a diagonal element portion of x.
[0080]
In step S118, it is determined whether i2 <1. If the determination in step S118 is no, in step S119, the inner product of the householder vector × α is stored in the upper triangular part of w2, and i1 = i2-1 is set. In step S120, it is determined whether i1 <1. If NO, w2 (i1, i2) = w2 (i1, i1) × (a (is + i2: n, is + i2) is determined in step S121. -1)^tXa (is + i2: n, is + i1-1)), i1 = i1-1 is calculated, and the process returns to step S120. If the determination in step S120 is yes, i2 = i2-1 is set in step S122, and the process returns to step S118. If the determination in step S118 is YES, i1 = blk-2 is set in step S123, and the expansion coefficient is calculated in a double loop. The upper side of the triangular matrix is determined from right to left, and calculation is performed so that they are stacked. This is equivalent to determining the coefficient by adding the expansion when the householder transformation is applied from the left. In step S124, it is determined whether i1 <1. If the determination in step S124 is no, i2 = blk is set in step S125, and it is determined in step S126 whether i2 <i1 + 1. If the determination in step S126 is no, the element on the upper side is determined from left to right in step S127. Use the previous coefficient. That is, calculate w2 (i1, i2) = w2 (i1, i2) + w2 (i1, i1 + 1: i2-1) × w2 (i1 + 1: i2-1, i2), i2 = i2-1 Return to step S126. If the determination in step S126 is YES, i1 = i1-1 is set in step S128, and the process returns to step S124. If the determination in step S124 is YES, the process proceeds to step S129, i2 = blk is set, and in step S130, it is determined whether i2 <1. If the determination in step S130 is no, in step S131, a coefficient α that is insufficient in the following loop is applied. First, i1 = i2-1 is set, and in step S132, it is determined whether i1 <1. If the determination in step S132 is no, in step S133, w2 (i1, i2) = w2 (i1, i2) × w2 (i2, i2), i1 = i1-1 is calculated, and the process returns to step S132. If the determination in step S132 is yes, i2 = i2-1 is set in step S134, and the process returns to step S130. If the determination in step S130 is yes, in step S135, the BU^tIs calculated and stored in W. W (1: blk, 1: iwidth) = TRU (w2) × W (1: blk, 1: iwidth) is calculated. The upper triangular matrix of U, the lower rectangular part, and the BU stored in W^tUse (1 + UBU^t) × EV is calculated. That is, 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) is calculated, and the process returns to step S115.
[0081]
According to the present invention, it is possible to provide a high-performance and scalable eigenvalue / eigenvector parallel calculation method using a shared memory type scalar parallel computer.
According to the embodiment of the present invention, in particular, the calculation relating to eigenvector conversion is about 10 times faster with a scalar machine than the conventional method. The eigenvalues and eigenvectors of real symmetric matrices using these algorithms are obtained by using the same Sturm method and inverse iteration method to obtain the eigenvalues and eigenvectors of the tridiagonal matrix. 6.7 times faster with 7 CPUs. Other routines of SUN (Inferior function: eigenvalues and eigenvectors cannot be obtained selectively) 2.3 times faster than those using the divide & conquer method for calculating eigenvalues and eigenvectors of tridiagonal matrices .
[0082]
The eigenvalues and eigenvectors of Hermitian matrices using these algorithms are obtained by using the same Sturm method and inverse iteration method to calculate the eigenvalues and eigenvectors of the tridiagonal matrix. 4.8 times faster with 7 CPUs than functions. Other SUN routines (poor in function: eigenvalues cannot be obtained selectively) 3.8 times faster than the divide & conquer method for calculating eigenvalues and eigenvectors of tridiagonal matrices.
[0083]
Please refer to the following textbooks for the basic algorithm of matrix calculation.
G. H. Golub and C. F. Van Loan "Matrix Computations" third edition The Johns Hopkins University Press, 1996
For parallel computation of tridiagonalization, refer to the following documents.
J. Choi, JJ Dongarra, and DW 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 Rid Laboratory managed by Martin Marietta Energy System, Inc., for the US DEPARTMENT OF ENERGY under Contract No. DE-AC05-84OR21400, ORNL / TM-12472. (Appendix 1) Parallel processing method of eigenvalue problem for shared memory type scalar parallel computer Is a program that realizes
Blocking a real symmetric matrix or Hermitian matrix to be computed, copying the blocked part to a memory work area, and tridiagonalizing using the product between the blocks;
Obtaining eigenvalues and eigenvectors from the tridiagonalized matrix;
Transforming the eigenvector obtained from the tridiagonalized matrix by Householder transformation modified to be a parallel operation of a matrix operation of a predetermined block width to obtain an eigenvector of the original matrix;
A program for causing a shared memory scalar parallel computer to implement a parallel processing method characterized by comprising:
[0084]
(Supplementary note 2) The program according to supplementary note 1, wherein in the tridiagonalization step, the block portion is updated by a recursive program.
(Supplementary Note 3) In the tridiagonalization step, the block portion is divided and stored in small blocks so as not to straddle the pages of the cache memory, and processing is performed by performing parallel operations on each processor. The program according to attachment 1, wherein the program is executed.
[0085]
(Supplementary note 4) In the step of obtaining the eigenvector, the matrix for the householder transformation is generated simultaneously by each processor as an upper triangular matrix that is a small coefficient matrix that can be processed by each processor. The program according to 1.
[0086]
(Supplementary Note 5) In the step of obtaining the eigenvector, the eigenvector divides the second dimension of the stored two-dimensional array evenly according to the number of processors, and assigns each divided area to each processor. The program according to supplementary note 1, wherein the program is calculated by:
[0087]
(Appendix 6) A parallel processing method for an eigenvalue problem for a shared memory type scalar parallel computer,
Blocking a real symmetric matrix or Hermitian matrix to be computed, copying the blocked part to a memory work area, and tridiagonalizing using the product between the blocks;
Obtaining eigenvalues and eigenvectors from the tridiagonalized matrix;
Transforming the eigenvector obtained from the tridiagonalized matrix by Householder transformation modified to be a parallel operation of a matrix operation of a predetermined block width to obtain an eigenvector of the original matrix;
A parallel processing method comprising:
[0088]
(Supplementary note 7) The parallel processing method according to supplementary note 6, wherein in the tridiagonalization step, the block portion is updated by a recursive program.
(Supplementary note 8) In the tri-diagonalization step, the block portion is divided into small blocks so as not to cross the pages of the cache memory, and is processed by performing operations in parallel on each processor. The parallel processing method according to attachment 6, wherein the parallel processing method is performed.
[0089]
(Supplementary note 9) In the step of obtaining the eigenvector, the matrix for the householder transformation is generated simultaneously in each processor by an upper triangular matrix that is a small coefficient matrix that can be processed by each processor. 6. The parallel processing method according to 6.
[0090]
(Supplementary Note 10) In the step of obtaining the eigenvector, the eigenvector divides the second dimension of the stored two-dimensional array evenly according to the number of processors, and assigns each divided area to each processor. The parallel processing method according to appendix 6, wherein the parallel processing method is calculated as described above.
[0091]
【The invention's effect】
Efficient and scalable eigenvalue / eigenvector solution can be realized.
[Brief description of the drawings]
FIG. 1 is a diagram showing a hardware configuration of a shared memory type scalar computer that is premised on an embodiment of the present invention.
FIG. 2 is a diagram (part 1) illustrating an algorithm according to an embodiment of the present invention.
FIG. 3 is a diagram (part 2) for explaining the algorithm of the embodiment of the present invention;
FIG. 4 is a diagram (part 3) for explaining the algorithm of the embodiment of the present invention;
FIG. 5 is a diagram (part 4) for explaining the algorithm of the embodiment of the present invention;
FIG. 6 is a diagram (part 5) for explaining the algorithm of the embodiment of the present invention;
FIG. 7 is a diagram (part 6) for explaining the algorithm of the embodiment of the present invention;
FIG. 8 is a diagram (part 7) for explaining the algorithm of the embodiment of the present invention;
FIG. 9 is a diagram (No. 8) for explaining the algorithm of the embodiment of the present invention;
FIG. 10 is a diagram (No. 9) for explaining the algorithm of the embodiment of the present invention;
FIG. 11 is a diagram (No. 10) for explaining the algorithm of the embodiment of the present invention;
FIG. 12 is a pseudo code (No. 1) of a routine according to the embodiment of the present invention.
FIG. 13 is a pseudo code (No. 2) of a routine according to the embodiment of the present invention.
FIG. 14 is pseudo code (No. 3) of a routine according to the embodiment of the present invention.
FIG. 15 is a pseudo code (No. 4) of a routine according to the embodiment of the present invention.
FIG. 16 is a pseudo code (No. 5) of a routine according to the embodiment of the present invention.
FIG. 17 is a pseudo code (No. 6) of a routine according to the embodiment of the present invention;
FIG. 18 is a pseudo code (No. 7) of a routine according to the embodiment of the present invention.
FIG. 19 is a first diagram illustrating a pseudo-code process in a flowchart.
FIG. 20 is a diagram (part 2) illustrating a process of pseudo code in a flowchart.
FIG. 21 is a diagram (part 3) illustrating pseudo code processing in a flowchart;
FIG. 22 is a (fourth) diagram illustrating a pseudo-code process in a flowchart.
FIG. 23 is a (fifth) diagram illustrating a pseudo-code process in a flowchart.
FIG. 24 is a sixth diagram illustrating a pseudo-code process.
FIG. 25 is a flowchart (No. 7) illustrating pseudo code processing;
FIG. 26 is a diagram (No. 8) illustrating pseudo-code processing in a flowchart.
FIG. 27 is a ninth diagram illustrating a pseudo-code process.
FIG. 28 is a diagram (No. 10) illustrating a pseudo-code process in a flowchart.
FIG. 29 is a flowchart (No. 11) illustrating pseudo-code processing.
[Explanation of symbols]
10-1 to 10-n processor
11-1 to 11-n memory module
12 Interconnection network
13-1 to 13-n Secondary cache memory

Claims (4)

A program for realizing a parallel processing method of an eigenvalue problem for a shared memory scalar parallel computer,
A real triangular matrix or Hermitian matrix A (1: k, 1: k) (k is an integer) to be calculated is represented by a column block A ((i−1) * iblk + 2: k, (i−1) * Iblk + 1: i * iblk) = U, remaining block A ((i−1) * iblk + 2: k, i * iblk + 1: k) = A _i (where i is the number of repetitions and iblk is the block width) Copying the block portion U into the work area of the memory and tridiagonalizing using the product between the blocks;
Obtaining eigenvalues and eigenvectors from the tridiagonalized matrix;
In each processor of the shared memory type parallel computer, eigenvectors EV (1: n, 1+ (k−1) * nb: k * nb) obtained from the tridiagonalized matrix (determined eigenvectors as nev In order to perform parallel computation by calculating matrix operations of a predetermined block width nblk on these divided eigenvectors, respectively, and assuming that nb is the number to be processed equally by p processors.
Householder transformation 1 + α _k * u _k u _k ^t (hereinafter, referred to as tridiagonalization of real symmetric matrix or Hermitian matrix) in the lower triangular part of A (1: n, 1: n−2) of the square matrix A for Hermitian matrix t → H: H in the transposed complex conjugate) of alpha _k diagonal elements a (k, k), also the Householder vector _{u k} a _(k + 1: storing n, k) the Every other (the block width obtained by equally dividing 1 to n−2 is nblk) and A (1: n, n−2−m * nblk + 1: n−2− (m−1) according to the repetition number m starting from 1. ) * Nblk) The bundle of householder vectors stored below the diagonal of the part is U _m, and a matrix B of nblk × nblk is obtained by each processor, and A (1: n, n−2−m * nblk + 1) : Household stored in n-2- (m-1) * nblk) From the smaller the loaders vector, multiplied by the alpha _i in the inner product of the corresponding Householder vector u _i and u _j expressed by subscript indicated by numbers renumbered from 1 to nblk i, j (i <j ) Is stored in the upper triangular part B (i, j), α _j corresponding to the householder vector corresponding to the number j that has been reassigned is stored in the diagonal element B (i, j) of B, and these data are Further, by updating in the upper triangular part, the same calculation is performed in each processor and the same matrix B is formed,
Each processor uses the matrix B to transform it into a householder transformation in the form of a matrix operation of (1 + U _m BU _m ^t ), and the eigenvector EV (1: n, 1+ (k−1) * nb: k * nb) Performing a parallel operation in each processor in the form of a matrix operation from the left and calculating these eigenvectors by repeating these processes with the number of repetitions m incremented by +1,
A program for causing a shared memory scalar parallel computer to implement a parallel processing method characterized by comprising:   The program according to claim 1, wherein in the tridiagonalization step, the block portion is updated by a recursive program.   In the tri-diagonalization step, the block portion is divided and stored in small blocks so as not to cross the pages of the cache memory, and each processor performs processing in parallel. The program according to claim 1, characterized in that:   In the step of obtaining the eigenvector, the eigenvector is calculated by equally dividing the second dimension of the stored two-dimensional array according to the number of processors and assigning each divided area to each processor. The program according to claim 1, wherein:

Images