JP5110081B2 - 共有メモリ型スカラ並列計算機向け、実対称行列の三重対角化の並列処理方法 - Google Patents

Description

本発明は、共有メモリ型スカラ並列計算機における、実対称行列の三重対角化を高速に処理可能な並列処理方法に関する。

実対称行列の固有地問題を解く上で対角化を行い、三重対角行列の固有値問題に変形して解く方法がある。この方法では、実対称行列を三重対角化する部分の計算コストが大きい。この三重対角化に関して、共有メモリ型の並列計算機では、ベクトル計算機のような強力なメモリアクセス能力がないため、メモリアクセスに対して演算量を増やすアルゴリズム上の工夫がなされている。

特許文献１には、ブロック化された三重対角化の処理の基本が記述されている。
しかし、計算のコストの大部分が行列ベクトル積であるため、メモリアクセスのスピードの影響を大きく受ける。全体の性能を改善するために、この部分の改良が必要である。

演算量を増やす工夫を行ったブロック化した方法では、行列ベクトル積を連続アクセスを使って計算すること、および、その計算を各ＣＰＵで均等に並列処理することがポイントになる。このために、更新時に必要となる行列全体を更新するとき、下三角部分を更新し、上三角部分をコピーしていた。更新時の負荷をバランスさせるために、行列を列方向にCPU の総数の 2倍で分割し、更新での演算量がバランスするようにi 番目と2 ×＃ＣＰＵ−（ｉ−１）番目をペアとして、各ＣＰＵに割り付けていた。しかし、この方法は、行列ベクトル積の計算で参照するメモリ領域が離れており、このためキャッシュ上のデータと干渉して、データをキャッシュに保存しにくい欠点があった。

行列ベクトル積、行列積による更新および負荷分散がうまく調和し、行列ベクトル積を高速に行なえる方法を探せれば、大きな性能向上を得られる可能性がある。
なお、非特許文献１には、三重対角化の基本アルゴリズムが、非特許文献２には、三重対角化の並列計算が記載されている。
特開2004-5528号公報 G.H. Golub, C.F. van Loan, Matrix Computation Second Edition, Johns Hopkins University Press 1989 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 U.S. DEPARTMENT OF ENERGY under Contract No. DE-AC05-84OR21400, ORNL/TM-12472.

本発明の課題は、共有メモリ型スカラ並列計算機において、高速に実対称行列の三重対角化演算を行うことが出来る並列処理方法を提供することである。
本発明の並列処理方法は、複数のプロセッサを備える共有メモリ型スカラ並列計算機において、実対称行列の三重対角化を高速に行うための並列処理方法であって、該実対称行列の下三角行列を、プロセッサの個数個の列ブロック行列に分割し、該列ブロックを各プロセッサの記憶領域にロードし、各プロセッサは、１回のロードで読み込んだ列ブロックの要素に対し、該実対称行列の縦方向についての演算と、横方向についての演算とを施すことを特徴とする。

本発明が念頭に置く共有メモリ型並列計算機のハードウェア構成を示すブロック図である。 ブロック化された三重対角化演算について説明する図（その１）である。 ブロック化された三重対角化演算について説明する図（その２）である。 ブロック化された三重対角化演算について説明する図（その３）である。 ブロック化された三重対角化演算について説明する図（その４）である。 ブロック化された三重対角化演算について説明する図（その５）である。 本発明の実施形態に従った並列処理方法を説明する図（その１）である。 本発明の実施形態に従った並列処理方法を説明する図（その２）である。 本発明の実施形態に従った並列処理方法を説明する図（その３）である。 本発明の実施形態の並列処理方法をより詳しく説明する図（その１）である。 本発明の実施形態の並列処理方法をより詳しく説明する図（その２）である。 本発明の実施形態の並列処理方法をより詳しく説明する図（その３）である。 本発明の実施形態の並列処理方法をより詳しく説明する図（その４）である。 本発明の実施形態に従った処理のフローである。 本発明の実施形態に従った処理のフローである。 本発明の実施形態に従った処理のフローである。 本発明の実施形態に従った処理のフローである。 本発明の実施形態に従った処理のフローである。 本発明の実施形態に従った処理のフローである。 本発明の実施形態に従った処理のフローである。 本発明の実施形態に従った処理のフローである。 本発明の実施形態に従った処理のフローである。

本発明では、特許文献１に記載の三重対角化の方法を基本に、これに改良を加えたものである。
本発明の実施形態においては、三重対角化を行なうとき、更新部分の演算密度を高めるために、ブロック化されたアルゴリズムを使う。行列ベクトル積部分の対称性を利用して、下三角部分のみを利用して行列ベクトル積を並列に実行する。このとき２回の行列ベクトル積が出現するが、メモリアクセスを共通にして１回の参照にする。下三角部分のみ使うため、更新を行なうための下三角行列を演算量が均等になるように分割を動的に設定する。これにより余分なコピーを行なう必要がなくなる。また、行列ベクトル積で参照したメモリがキャッシュに保持されやすくなる。

下三角部分の行列積(rank-p update) による更新の計算量だけでなく、行列を参照する行列ベクトル積の演算量も同じく均等( 下三角行列を分割したおのおのの面積が均等）になるように列方向に分割して処理する。このとき、分割された領域のおのおのは下三角形部分と長方形部分の 2つの領域からなると見なせる。下三角形部分の更新は、この部分を適当に再分割し、いくつかの行列積で部分的に更新する。対角要素を含む分割された下三角部分を通常の行列積で更新すると、上三角部分の一部も更新するが、なるべくはみだす部分の計算を少なくするようにすることで、演算量を削減し、下三角部分だけを更新するときの計算に近い演算量で更新を行なう。

図１は、本発明が念頭に置く共有メモリ型並列計算機のハードウェア構成を示すブロック図である。
図１に示されるように、本発明が念頭に置く並列計算機では、ＣＰＵ１２−１〜１２−ｍとメモリモジュール１０−１〜１０−ｎがバスで構成される相互結合網１１によって相互に接続された構成となっている。ＣＰＵ１２−１〜１２−ｍは、それぞれ、L2キャッシュ及びバスインタフェースと、これに接続される複数のL1キャッシュとＣＰＵコアからなっている。この構成をマルチコアＣＰＵという。マルチコアＣＰＵは、１つのＣＰＵの内部にＣＰＵコア（L1キャッシュも内臓）を複数封入したもので、各コアからL2キャッシュ及びバスインタフェースは共通に使う構成である。論理的には、各ＣＰＵコアが１つのＣＰＵとして動作するように見える。図１の例では、１つのＣＰＵに２つのＣＰＵコアが封入された例を示しているが、４つ封入されるものもある。

図２〜図４Ｃは、ブロック化された三重対角化演算について説明する図である。
図２〜図４Ｃを参照して、固有値・固有ベクトルを求めるための三重対角化の数学低枠組みを説明する。
ａ）ブロック化された三重対角化の数学的アルゴリズム
行列をブロック幅ごとに３重対角化する。行列をブロックに分割してブロック単位で下記のブロックアルゴリズムで三重対角化を行う。図２はn番目のブロックを処理しているときの図である。

最後のブロックは、そのブロック幅−２をブロックとする左詰めのブロックに関してアルゴリズムを適用してすべての処理を終わる。
step1:A_n のn+i 行ベクトルよりHouseholder ベクトルuを作る。
step2: v_i = A_{n +i}u , w_i = v_i - u(u^tv)/2, を計算する。
step3: U_i =(U_i-1,u_i ) , W_i =(W_i-1,w_i )と更新する。
step4: if(i < blks)then
A_nのn+i+1 列目を更新する。

A_n(*,n+i+1) = A_n (*,n+i+1)- U_i W_i (n+i+1,*)^t- W_i U_i (n+i+1,*)^t
endif
enddo
step５：A_n+blks = A_n - U_blks W_blks ^t- W_blks U_blks ^t
次に、ブロック化されたHouseholder 変換による三重対角化について説明する。ブロック化されたHouseholder変換は、以下の式の通りに行う。
ｖ＝（v₁,v₂,...,v_n )
｜ｖ｜² ＝ｖ・ｖ＝ｈ² 、
U_v =(h,0,…,0)^t とすると、U_v = v − （v₁−h ，v₂，．．．，v_n）の関係がある。
U ＝（I −uu^t / ｜u ｜² ）＝（I −αuu^t ）, ここでu=(v₁−h ，v₂，．．．，v_n )
以下αを省略して計算する。
A_n+1 ＝U^t A_n U ＝（I-uu^t ）A （I-uu^t）
＝A_n - uu^t A_n - A_n uu^t＋ uu^t A_n uu^t
＝A_n - uv^t -vu^t ＋uv^t uu^t
＝A_n - uw^t -uu^t u^t v/2 - wu^t- uu^t u^t v/2 +uu^t u^t v
＝A_n - u w^t -wu^t ・・・（＊）
ここで、 w = v - u(u^t v)/2, v = A_n u
これを繰りし利用すれば、
A_n+k ＝A_n −U_k W_k ^t −W_kU_k ^t ・・・（＊＊）
ｋステップ目の計算は（＊）および（＊＊）からV_n を次のように計算できる。
v_k = A_n u^k -U_k-1 W_k-1 ^tu_k - W_k-1U_k-1 ^t u_k
w_k =v_k - u_k u_k ^t v_k/2
U_k =(U_k-1,u_k ) , W_k =(W_k-1,w_k)
A_n+k = A_n - U_k W_k ^t - W_kU_k ^t
ｂ）Householder変換を構成する情報の格納
固有ベクトルを計算するときに３重対角化でつかったHouseholder 変換が必要になる。このため、Ｕn およびαをhouseholder 変換を構成するベクトルの位置に格納する。αは対応する対角要素の位置に格納する。
ｃ）Ｕi を効率的に求める方法
ブロック部分の三重対角化を行うためのHouseholder 変換のために、次のベクトルを更新する必要がある。これらの計算をできるだけローカルに行うために、ブロック幅部分を作業域にコピーして３重対角化を行い、もとの領域に格納する。次の列ベクトルの更新をそのたび行わず、行列積の形にして演算密度を上げて計算を行う。このため、ブロック部分の三重対角化を再帰的なプログラミングで行う。

すなわち、図３の左の図における行列Ａのブロック領域を、図３右の作業領域Ｕにコピーし、当該ブロック領域を前半更新部分と後半更新部分とに分けて、再帰的に処理する。図４Ａ〜図４Ｃは、再帰的な処理の様子を示している。今、図４Ａの斜線部分を更新したとすると、次に、図４Ｂの斜線部分を更新し、更に、図４Ｃの斜線部分を更新する。このように、前半部分と後半部分とに分けて更新する仕方を入れ子状に適用する。すなわち、再帰的プログラムが深さ２まで呼ばれると、最初の前半処理で図４Ａの斜線部分をＢとして更新し、次に、図４Ｂの斜線部分、最後に図４Ｃの斜線部分が更新される。更新時の並列化は更新部分ブロック行列を行ベクトル方向に均等に分割して、各部分の更新をプロセッサで均等に分割した領域に対して並列して行う。

図５〜図７は、本発明の実施形態に従った並列処理方法を説明する図である。
１）行列ベクトル積を下三角部分のみを使って１つのアクセス（ロード）で 2倍の演算を行ない、各ＣＰＵで均等に行なう方法
図５に示されるように、下三角行列部分のみを利用して行列ベクトル積を計算する。行列ベクトル積で使う行列の大きさをＮ_ｎ×Ｎ_ｎとする。４並列で行うとき、下三角部分を４つの部分に分割する位置は、以下のようになる。
i_m=N_n×(1-sqrt(1-m/#p))
ｍ＝１、２、３で計算したものがちょうど境界位置になる。なお、＃ｐは、分割する総数、つまり、並列数である。

Ａ_nの対称性つまり、Ａ_n ＝ Ａ_n ^t を利用して、各プロセッサで自分の受け持つ領域にある情報から、行列積を計算する。
各部分は対角ブロック行列部分の下三角行列部分Ｄ_m と、その下にある長方形の行列Ｒ_mとからなる。

下三角行列を#Pに分割した、m 番目の部分を転置したものは、ちょうど上三角部分を行方向に分割したものになっている。この関係から、以下の計算を行なうことができる。
ｖ_xm1 (i_m＋１：Ｎ_n )=Ｒ_m ×ｕ₁(ｉ_m-1＋１：ｉ_m) （式１）
ｖ_xm2 (i_m-1＋１：ｉ_m )=Ｒ_m ^t ×ｕ₁(ｉ_m＋１：Ｎ_n) （式２）
ｖ_xm3 (i_m-1＋１：ｉ_m )=Ｄ_m ×ｕ₁(ｉ_m-1＋１：ｉ_m) （式3 ）
ｖ_xm4 (i_m-1＋１：ｉ_m )=ｎｏｔｄｉａｇ（Ｄ_m ）^t×ｕ₁(ｉ_m-1＋２：ｉ_m ) （式4 ）
（ｎｏｔｄｉａｇ（Ｘ）は、下三角行列Ｘの対角部分を除いた下三角行列を表す）

この計算で、行列の要素を連続にアクセスしながら、１つの要素のアクセスに対して、ｖ_xm1とｖ_xm2の計算を同時に行なうことで、行列ベクトル積の性能を 2倍にすることができる。

式１と式２で同時に計算するプログラム例は以下のようになる。
簡単のためＲ_mを行列ｎ×ｍの行列、ｘ_１を長さｎのベクトル、ｘ_２を長さｍのベクトル、ｙ_１を長さｍのベクトル、ｙ_２を長さｎのベクトルとする。

ｘ_１＝Ｒ×ｙ_１、および、ｘ_２＝Ｒ^t ×ｙ_２ をＲの参照を共通に行なう。

x1(1:n)=0
do j=1,m
tmp1=y1(j)
sum=0
do i=1,n
tmp2=y2(i)
tmp=R(i,j)
x1(i)=x1(i)+tmp*tmp1
sum=sum+tmp*tmp2
enddo
x2(j)=sum
enddo
この場合、x1やy2はデータ量が小さく、かつ、繰り返しアクセスされるので、キャッシュに保持され、高速にアクセスできる。R(i,j)のアクセスに比べ、演算はほとんど無視できるほど速いので、別々に行列ベクトル積を行なう場合に比べて約２倍の性能になる。

（＊）分割の境界を計算する式、i_m = N_n ×(1-qsrt(1-m/#p))は次のように求めることができる。
各CPU での演算数が等しい。つまり、要素数に比例する行列の対応部分の面積（あるいは要素数そのもの）が等しいと考える。
Ｎ_n ²-( Ｎ_n - ｉ_m )²＝Ｎ_n ²×m/#p，ｒ＝ｉ_m ／Ｎ_n と置いて計算する。

ｒ² −２ｒ＝m/#p ， この根の0<= ｒ<=1 をとる。
(ｒ−１)²＝m/#p
ｒ＝１−sqrt(1-m/#p) ( 根の範囲より小さい方を採用する。）
２）更新部分の並列化および対角ブロック部分の計算方法
step5の更新部分は対称性を利用して下半分を計算する。更新の演算量が各CPU(スレッド) で、均等化するように、分割する。分割方法は行列ベクトル積で分割したのと同じ方法で行なう。

行列ベクトル積で使う行列の大きさをN_n ×N_n とする。4 並列で行なうとき、下三角部分を 4つの部分に分割する位置は以下のようになる。
i m = N_n ×(1-qsrt(1-m/#p)) m=1,2,3 で計算したものがちょうど境界の位置になる。#pは、分割する総数、つま並列数である。

このように、ＣＰＵごとに分割された領域は、対角ブロック行列部分とその下の長方行列部分に分かれる。図６は、配列A_n+iをＣＰＵごとに分割した様子を示している。長方行列部分は、普通の行列積で
A_n+k ＝A_n −U_k W_k ^t −W_kU_k ^t から、対応部分を計算する。
対角ブロック行列部分の更新は、図７のように行なう。すなわち、対角ブロックの下三角部分について、対角線の２等分点を決める。この等分点から、行方向、及び、列方向に下三角行列を分割する。矩形部分は、行列積で更新する。２つの大きさが半分の下三角行列ができるので、これらに関して、同じく対角線上の等分点を決め、同様に列方向、行方向に分割する。そして、対象となる下三角行列の大きさが十分小さくなったら、上三角行列も一緒に行列積で計算する。これらは、簡単な再帰プログラムで作ることが出来る。

以上の並列処理方法により、高性能かつスケーラビリティのある実対称行列の三重対角化を実現できる。以上の方法は、従来の方法に比べて２倍超の性能であることがわかった。

図８〜図１１は、本発明の実施形態の並列処理方法をより詳しく説明する図である。
ブロック化された三重対角化演算においては、以下のような処理を行う。
１）三重対角化でブロックごとに、三重対角化を行う。上のブロックを三重対角化し、その三重対角化で生成された行列U およびW を使って、ブロック幅分縮小された正方行列( この下三角行列は、図８において、太い点線で表示) を更新する。ただし、この正方行列は対称行列なので、下三角行列部分( 太い点線で囲んだ三角形部分) のみ更新して上三角行列部分のデータとして利用する。
２）ブロックごとの三重対角化では、ブロック内の太線で示した正方行列(A_n+1：太線で囲んだ正方行列) に関する行列ベクトル積を計算する必要がある。この計算で、正方行列の大きさは１ずつ小さくなっていく。この部分も対称行列なので、下三角行列部分を使って計算する。

これらの計算が計算の中での演算コストが大きい。並列処理で効果を引き出すには以下の点の考慮が必要である。
・メモリアクセスが遅いスカラ計算機で行列ベクトル積の計算を行うとき、ほとんどがメモリアクセスの時間である。このた、メモリアクセスのコストを削減することがポイントである。
・行列の対称性を利用して、下三角行列部分を使うとき、並列処理の観点から負荷分散が均等であることが望ましい。

そこで、以下のような構成を有するプログラムを用意する。
図９は、本発明の実施形態に従ったプログラムの機能ブロック図である。
ブロック三重対角制御部１５は、ブロック部分の三重対角化部１６と、更新部１９からなる。三重対角化部１６は、各スレッドが担当する下三角行列の２次元目の区間を決める行列ベクトル積の負荷制御部１７を備え、負荷制御部１７は、ブロック対角行列の下三角部分とその下の長方形行列で計算を行う行列ベクトル積並列計算部１８を備える。また、更新部１９は、各スレッドが担当する下三角行列の２次元目の区間を決める行列積による演算の負荷制御部２０を備える。負荷制御部２０は、更に、ブロック対角部分を除く更新を行うブロック更新部２１とブロックの細分化の制御を行うブロック対角部の下三角行列部の更新制御部２２からなり、更新制御部２２は、ブロック対角行列を細分化して下三角行列を細かい正方行列で近似して計算する計算部２３を備える。

図１０は、下三角行列（繰り返しごとに大きさが１ずつ小さくなる）から、行列ベクトル積を計算する処理を説明する図である。

ここでは、V(nptr:n) = A(nptr:n, nptr:n) * U(nptr:n)の計算をする。ここで、nptrは、下三角行列の左端、nは、行列の右端、msは、各スレッドの受け持つ始点、meは、各スレッドの受け持つ終点を示す値である。

配列として、V(*,i)、U(*,i)、vsum(1:n,numthrd)と、shared属性の全体行列Aを確保する。そして、以下の処理を行う。
１）vsum(nptr:n,nothrd)=0と各スレッドで０クリアする。
２）以下の２つの計算をL2の列方向へ、図１０のL2の中の矢印のように参照しながら同時に計算する。後に説明するmatvec1のフローを参照されたい。
vsum(me+1:n)=L2*U(ms:me)
vsum(ms:me)=L2^t*U(me+1:n)
３）L1を使ってブロック対角部分を同様に計算する。
vsum(ms:me)=vsum(ms:me)+L1*U(ms:me)+ND(L1)^t*U(ms:me,i)
（ここで、ND(L1)は、L1の対角要素に０とした下三角行列）
４）nptr:nをスレッドごとに均等分割した区間を作る。isを始点、ieを終点とする。計算終了をバリア同期をとって確認し、各スレッドで並列計算する。
V(is:ie,i)=Σvsum(is:ie,nothrd)
nothrd=1〜numthrd
（nothrdは、スレッド番号１〜numthrdの値を取る。numthrdは、総スレッド数。）
図１１は、更新部の処理を説明する図である。

A(nptr:n,nptr:n)の下三角行列部分を以下の式で更新する。
A(nptr:n,nptr:n)=A(nptr:n,nptr:n)-U(nptr:n,1:blk)*W(nptr:n,1:blk)^t-W(nptr:n,1:blk)*U(nptr:n,1:blk)^t
ここで、nptrは、下三角行列の左端、blkは、ブロック幅、nは、行列の右端、msは、各スレッドの受け持つ始点、meは、各スレッドの受け持つ終点である。

L1を、底辺と垂線の中点を求めて、L11、L12、L22に分割する。更に、L11、L22を同様に、中点を求めて分割し、下三角行列を正方行列の集合に分割する。これを繰り返す。ブロック対角行列が十分小さくなったら、一般の行列の行列積で計算する。これは、再帰プログラムで実現する。（詳細は、サブルーチンltgemmtrtのフローを参照。）
図１２〜図２０は、本発明の実施形態に従った処理のフローである。

図１２は、実対称行列を三重対角化するサブルーチンtridのフローである。このサブルーチンでは、shared配列A(k,n)、diad(n)、及び、sdiag(n)を入力とする。diag、sdiagは、計算した三重対角行列の対角要素、副対角要素を出力として返却するための配列である。

ステップＳ１０において、作業域U(n+1,iblk)、v(n+1,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を設定する。ステップＳ１４において、作業領域Uにブロック三重対角化の対象領域をコピーする。すなわち、U（nbase+1:n,1:iblk)←A(nbase+1:n,nbase+1:nbase+iblk)とする。ステップＳ１５において、サブルーチンblktridを呼び出し、Uにコピーした部分の三重対角化を行う（istart=1、ブロック幅はiblkを受け渡す）。ステップＳ１６において、三重対角化されたものを配列Aに戻す。すなわち、A(nbase+1:n,nbase+1:nbase+iblk)←U(nbase+1:n,1:iblk)とする。ステップＳ１７において、サブルーチンupdateを呼び出し、A(nbase+iblk:n,nbase+iblk:n)の下三角部分を更新する。そして、ステップＳ１２に戻る。

ステップＳ１８では、nbase=(nb-1)*iblk、istart=1、iblk2=n-nbaseと設定する。ステップＳ１９において、作業領域Uにブロック三重対角化の対象領域をコピーする。すなわち、U(nbase+1:n,1:nwidth)←A(nbase+1:n,nbase+1:n)とする。ステップＳ２０において、サブルーチンblktridを呼び出し、Uにコピーした部分の三重対角化を行う（istart=1、ブロック幅はiblk2を受け渡す）。ステップＳ２１において、三重対角化されたものを配列Aに戻す。すなわち、A(nbase+1:n,nbase+1:n)←U(nbase+1:n,1:width)とする。ステップＳ２２において、並列処理のために生成したスレッドを消し、このサブルーチンを抜ける。

図１３は、サブルーチンblktridのフローである。
このサブルーチンは、再帰的プログラムである。呼び出すプログラム文は、subroutine blktrid ( A ,k ,n ,diag ,sdiag ,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, ぇん(ｎ-iptr+numthrd-1)/numthrd,
is=iptr+(nothrd-1)*len+1、 ie=min(n,iptr+nothrd*len)
ステップＳ３１において、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)*U(is3:ie3,is2:ie2)^tとする。

ステップＳ３２において、スレッド間でバリア同期を取り、ステップＳ３３において、サブルーチンblktridを再帰的に呼び出して、このサブルーチンを抜ける。
図１４は、サブルーチンbtunitのフローである。

このサブルーチンの呼び出しは、btunit(A, k, n, diag, sdiag, nbase, istart, nwidth, U, V, nothrd, numthrd)とする。なお、以下で、tmpという配列は、計算用一時配列である。ステップＳ３５において、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)=U(is:ie,i)^t*U(is:ie,i)とし、ステップＳ４２で、バリア同期を取る。ステップＳ４３において、nothrd=1であるか否かを判断する。ステップＳ４３の判断がＮｏのときは、ステップＳ４５に進む。ステップＳ４３の判断がＹｅｓの場合には、ステップＳ４４において、各スレッドで部分計算したデータの和の平方根をとり、三重対角化のための計算を行う（ハウスホルダーベクトルの作成）。すなわち、以下の計算を行う。なお、sumは、和をあらわし、sqrtは、平方根をあらわす。
sigma=sqrt(sum(tmp(1:numthrd)))
diag(iptr2)=u(iptr2,i), sdiag(iptr2)=-sigma,
U(nbase+i+1,i)=U(nbase+1,i)+sign(u(nbase+i+1,i)*sigma
alpha=1.0/(sigma*u(nbase+i+1,i), U(iptr2,i)=alpha
ステップＳ４５では、バリア同期を取る。ステップＳ４６において、iptr3=iptr2+1とする。ステップＳ４７において、A, V, U, n, k, I, iptr3, is, ie, nothrd, numthrdを受け渡して、サブルーチンmatvecを呼び出す。すなわち、A(iptr3:n,iptr3:n)の下三角行列からV=A*Uの計算を行う。ステップＳ４８において、バリア同期を取る。ステップＳ４９において、以下の計算をする。
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)))
ステップＳ５０において、バリア同期を取る。ステップＳ５１において、tmp(nothrd)=V(is:ie,i)^t*U(is:ie,i)とし、ステップＳ５２において、バリア同期を取る。ステップＳ５３において、nothrd=1であるか否かを判断する。ステップＳ５３の判断がＮｏの場合には、ステップＳ５５に進む。ステップＳ５３の判断がＹｅｓの場合には、ステップＳ５４において、beta=0.5*alpha*sum(tmp(1:numthrd))とする。ここで、sumは、ベクトルの和である。ステップＳ５５で、バリア同期を取り、ステップＳ５６において、V(is:ie,i)=V(is:ie,i)-beta*U(is:ie,i)とし、ステップＳ５７において、バリア同期を取る。ステップＳ５９において、ptr2<n-2であるか否かを判断する。ステップＳ５９の判断がＹｅｓの場合には、ステップＳ６０において、U(is:ie, i+1)=U(is:ie, i+1)- U(is:ie, istart:i)*V(i+1, istart:i)^t- V(is:ie, istart:i)*U(n+1, istart:i)^tを計算し、ステップＳ３８に戻る。ステップＳ５９の判断がＮｏの場合には、ステップＳ６１において、、U(is:ie, i+1:i+2)=U(is:ie, i+1:i+2)- U(is:ie, istart:i)*V(i+1:n, istart:i)^t- V(is:ie, istart:i)*U(n+1:n, istart:i)^tを計算し、サブルーチンを抜ける。

図１６は、サブルーチンmatvecのフローである。
引数は、A, V, U, n, k, I, is, ie, iptr3, nothrd, numthrdである。ステップＳ６５において、vsum(n,numthrd)をshared属性で確保する。ステップＳ６６において、各スレッドでvsum(iptr３:n,nothrd)=0と初期化する。制御部では、行列ベクトル積を各スレッドで分担計算するために、下三角行列A(iptr3:n, iptr3:n)を各スレッドで分担する区間を計算する。すなわち、以下の計算をする。
nn=n-iptr3+1,
xleft=1.0-dble(nothrd-1)/dble(numthrd),
xright=1.0-dble(nothrd)/dble(numthrd),
ms=nbase2+nn*(1.0-dqsrt(xleft))+1,
me=nbase2+nn*(1.0-dsqrt(xright))
ステップＳ６８において、matvec1を呼び出す。下三角行列の２次元目の範囲が始点ms、終点meの対角行列部分の下のA(me+1:n, ms:me)から各スレッドで
vsum(ms:me, nothrd)=A(me+1:n, ms:me)^t*U(me+1:n,i)
vsum(ne+1:n,nothrd)=A(me+1:n,ms:me)*U(ms:me)
を計算する。引数は、A, vsum, U, I, ms, me, is, ie, nothrd, numthrdである。

ステップＳ６９において、matvec2を呼び出す。下三角行列の２次元目の範囲が始点ms、終点meの対角行列部分の下三角行列から各スレッドで
vsum(ms:me)=vsum(ms:me)+LOW(A(ms:me, ms:me))*U(ms:me,i)
+NDLOW(A(ms:me, ms:me))^t*U(ms:me,i)
を計算する。

LOWは、対角要素を含む下三角行列部分であり、NDLOWは、対角要素を０とした下三角行列部分である。引数は、A, vsum, U, I, ms, me, is, ie, nothrd, numthrdである。
ステップＳ７０において、バリア同期を取り、ステップＳ７１において、

を計算し、このサブルーチンを抜ける。
図１７は、サブルーチンmatvec1のフローである。

ステップＳ７５において、k2=msとする。ステップＳ７６において、k2<=meであるか否かを判断する。ステップＳ７６の判断がＮｏの場合には、このサブルーチンを抜ける。ステップＳ７６の判断がＹｅｓの場合には、ステップＳ７７において、k1=me+1, sum=0.0, tmp1=U(k2,i)とする。ステップＳ７８において、k1<=nであるか否かを判断する。ステップＳ７８の判断がＹｅｓの場合には、ステップＳ７９において、tmp2=U(k1,i)、tmp=A(k1,k2)（Aの要素の１回のロードで２回演算する）、vsum(k1)=vsum(k1)+tmp*tmp1、sum=sum+tmp*tmp2、k1=k1+1とし、ステップＳ７８に戻る。ステップＳ７８の判断がＮｏの場合には、ステップＳ８０において、vsum(k2)=sum、k2=k2+1とし、ステップＳ７６に戻る。

図１８は、サブルーチンmatvec2のフローである。
ステップＳ８５において、k2=msとする。ステップＳ８６において、k2<=meであるか否かを判断する。ステップＳ８６の判断がＮｏの場合には、このサブルーチンを抜ける。ステップＳ８６の判断がＹｅｓの場合には、ステップＳ８７において、k1=k2+1、tmp1=U(k2,i)、tmpz=A(k2,k2)、sum=tmp1*tmpzとする。ステップＳ８８において、k1<=meであるか否かを判断する。ステップＳ８８の判断がＹｅｓの場合には、ステップＳ８９において、tmp2=U(k1,i)、tmp=A(k1,k2)（１回のロードで２回計算する）、vsum(k1)=vsum(k1)+tmp*tmp1、sum=sum+tmp*tmp2、k1=k1+1として、ステップＳ８８に戻る。ステップＳ８８の判断がＮｏの場合には、ステップＳ９０において、vsum(k2)=vsum(k2)+sum、k2=k2+1として、ステップＳ８６に戻る。

図１９は、サブルーチンupdateのフローである。
ステップＳ９５において、スレッド間でバリア同期を取る。ステップＳ９６は負荷の制御部で行う。ステップＳ９６においては、各スレッドでペアをつくり、更新を分担する始点、終点を決める。すなわち、以下の演算を行う。
nbase2=nbase+iblk,
nn=n-nbase2,
xleft=1.0-dble(nothrd-1)/dble(numthrd),
xright=1.0-dble(nothrd)/dble(numthrd),
is=nbase2+nn*(1.0-dqsrt(xleft))+1,
ie=nbase2+nn*(1.0-dsqrt(xright))
ステップＳ９７においては、以下を計算する。
A(ie+1:n, is;ie1)=A(ie+1:n, is:ie)-W(ie+1:n, 1:blk)*U(is:ie, 1:blk)^t-U(ie+1:n, 1:blk)*W(is:ie, 1:blk)^t
ステップＳ９８において、下三角行列部の更新（A=A-W*U^t)を行う。サブルーチンltgemmtrtを呼び出し、対角ブロック行列の下三角行列部分を更新する。引数として、W(is:ie, 1:blk), U(is:ie, 1:blk), A(is:ie, is:ie)及び対角ブロック行列の大きさlen=ie-is+1を渡す。

ステップＳ９９において、下三角行列部の更新(A=A-U*W^t)を行う。サブルーチンltgemmtrtを呼び出し、対角ブロック行列の下三角行列部分を更新する。引数として、U(is:ie, 1:blk), W(is:ie, 1:blk), A(is:ie, is:ie)及び対角ブロック行列の大きさlen=ie-is+1を渡す。すなわち、UとWを入れ替えて呼び出す。

ステップＳ１００において、スレッド間でバリア同期を取り、このサブルーチンを抜ける。
図２０は、サブルーチンltgemmtrtのフローである。

このサブルーチンは、対角行列の下三角行列部分の更新を行うもので、再帰的プログラムとして呼び出される。引数としては、W(is:ie, 1:blk), U(is:ie, 1:blk), A(is:ie, is:ie)をW(len,blk), U(len, blk), A(len, len)として受け取る。lenは、行列の大きさである。ブロックをどんどん小さく分割していく際にブロックが所定の最小値まで小さくなったら、ブロックをそれ以上小さく分割する処理をやめる。たとえば、ここでは、ブロックの最小値をnwc=8と設定しておく。ステップＳ１０５において、len<=nwcであるか否かを判断する。ステップＳ１０５の判断がＹｅｓの場合には、ステップＳ１０６において、A(len, len)=A(len, len)-W(len, blk)*U(len, blk)^tを計算し、このサブルーチンを抜ける。ステップＳ１０５の判断がＮｏの場合には、ステップＳ１０７において、len1=len/2、len2=len-len1、nptr=len1+1とし、ステップＳ１０８において、サブルーチンltgemmtrtを呼び出し、対角ブロック行列の下三角行列部分を更新する。このとき、引数として、U(1:len1, 1:blk), W(1:len, 1:blk), A(1:len, 1:len1)及び、対角ブロック行列の大きさlen1を渡す。

ステップＳ１０９において、A(nptr:len, 1:len2)=A(nptr:len, 1:len2)-W(nptr:len, blk)*U(1:len2, blk)^tを計算する。ステップＳ１１０において、len1=len2であるか否かを判断する。ステップＳ１１０の判断がＹｅｓの場合には、ステップＳ１１１において、len3=len2、nptr2=nptrとし、ステップＳ１１３に進む。ステップＳ１１０の判断がＮｏの場合には、ステップＳ１１２において、len3=len1、nptr2=nptr+1とし、ステップＳ１１３に進む。ステップＳ１１３では、サブルーチンltgemmtrtを呼び出し、対角ブロック行列の下三角行列部分を更新する。引数としては、U(nptr2:len, 1:blk), W(nptr2:len, 1:blk), A(nptr2:len, nptr2:len)と対角ブロックの大きさlen3を渡す。そして、ステップＳ１１３の後、このサブルーチンを抜ける。

Claims (2)

1. 共有メモリ型スカラ並列計算機において、ｎ×nなる実対称行列Aの三重対角化をブロック化した方法で行う方法であって、
   実対称行列を２次元配列に格納した場所に下三角行列部分だけを格納し、
   該実対称行列を幅iblkを持つ列ブロックごとに分割し、
   列ブロックごとに三重対角化を行うブロックの三重対角化部分と、三重対角化を行うときに生成されたブロック幅iblkの本数のベクトルを連結したブロック行列U,Wを使ってUW ^T とWU ^T を計算して三重対角化が完了したm番目の列ブロックまでを取り除いた残りのA(m*iblk+1:n,m*iblk+1:n)なる部分実対称行列の下三角行列を更新し、
   次の列ブロックの三重対角化を行うことを繰り返して、行列全体の三重対角化を行い、
   m番目の列ブロックに対する三重対角化部分で三重対角化を該列ブロックにあるi列に対して順次行うとき、i+1列を最左端列に持つ部分実対称正方行列A(i+1:n,i+1:n)とベクトルとの行列ベクトル積を行う場合に、A(i+1:n,i+1:n)の下三角行列を列方向に分割して各CPUに割り当てて該行列ベクトル積を並列に計算し、
   各CPUで、割り付けられた部分の行列の要素a(s,t)を順次取り出してこの値を使って、行列ベクトル積を行う場合、対角要素以外は対称なa(t,s)に対する行列ベクトル積の積を同じ値を使って、ベクトルの対応する要素と演算を行い、
   結果を最初にゼロクリアされた一次元配列の対応位置にそれぞれ加え、
   行列ベクトル積の各CPUでの並列計算の後、それぞれのCPUで計算した中間結果を加え合わせ、
   部分行列A(i+1:n,i+1:n)の下三角行列を列方向に分割して各CPUに割り当てるとき、A(i+1:n,i+1:n)が十分大きなとき、各CPUに割り当てられた下三角行列の要素数がおおよそ等しくなるように列方向に分割し、行列の列を各CPUに割り当てる、
   ことを特徴とする方法。
2. 列ブロックの三重対角化の処理の後ごとに行う更新部分に関しては、UW ^T ,WU ^T を計算して該更新部分から引くことで更新を行い、
   更新を行う下三角部分行列の更新に関して、この行列が十分大きなとき、各CPUに更新する要素数がおおよそ等しくなるようにこの下三角行列を列方向に分割して各CPUに割り当て、
   各CPUに割り当てられた、列を束ねたブロックが、対角ブロック部分が下三角行列と、その下の部分の長方形の部分に分かれる場合に、長方形の部分は行列積で計算し、
   対角ブロック部分の下三角行列部分は、対角ブロックの対角要素のほぼ中央で、対角ブロックを列方向、行方向に、もとの下三角行列の1/4の大きさの２つの下三角行列と長方形の行列部分に分割して３つの部分に分け、
   このように分割してできた下三角行列に対して、該分割を再帰的に適用して、生成された下三角行列が十分小さくなるまで繰り返して、長方形部分は行列積で更新し、対角要素を含む小さな下三角行列はこれと対称な位置関係にある上三角部分も合わせた正方行列を行列積で計算して更新する、
   ことを特徴とする請求項１に記載の方法。

Images