(57)【要約】  実対称行列の三重対角化を共有メモリ型スカラ並列計算機で行う場合、行列の更新を各スレッドに分割して行う。このとき、行列が対称行列であるので、下三角行列部分のデータを使って計算する。下三角行列部分は、縦に分割し、各ＣＰＵに割り当てるが、その際、各ＣＰＵに割り当てる要素数が同じになるように割り当てる。そして、１回のデータのロードで、縦方向の演算と横方向の演算を同時に行い、メモリロードの回数を減らす。対角ブロック行列部分の下三角行列部分は、１つの小さな正方行列と２つの小さな下三角行列に分割することを再帰的に行いながら更新をする。
2008136045.tif 000002
【特許請求の範囲】
【請求項１】
複数のプロセッサを備える共有メモリ型スカラ並列計算機において、実対称行列の三重対角化を高速に行うための並列処理方法であって、
該実対称行列の下三角行列を、プロセッサの個数個の列ブロック行列に分割し、該列ブロックを各プロセッサの記憶領域にロードし、
各プロセッサは、１回のロードで読み込んだ列ブロックの要素に対し、該実対称行列の縦方向についての演算と、横方向についての演算とを施す
ことを特徴とする並列処理方法。
【請求項２】
前記実対称行列の対角ブロック行列の下三角行列を更に小さい１つの正方行列と２つのした三角行列に分割する処理を再帰的に行いながら、該対角ブロック行列の更新を行うことを特徴とする請求項１に記載の並列処理方法。
【請求項３】
前記実対象行列の列ブロックへの分割は、各プロセッサに割り当てられる要素数が同じになるように行うことを特徴とする請求項１に記載の並列処理方法。
【請求項４】
前記プロセッサのそれぞれは、複数のＣＰＵコアを供えることを特徴とする請求項１に記載の並列処理方法。
【請求項５】
複数のプロセッサを備える共有メモリ型スカラ並列計算機において、実対称行列の三重対角化を高速に行うための並列処理装置であって、
該実対称行列の下三角行列を、プロセッサの個数個の列ブロック行列に分割し、該列ブロックを各プロセッサの記憶領域にロードするロード手段と、
各プロセッサは、１回のロードで読み込んだ列ブロックの要素に対し、該実対称行列の縦方向についての演算と、横方向についての演算とを施す演算手段と、
を備えることを特徴とする並列処理装置。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、共有メモリ型スカラ並列計算機における、実対称行列の三重対角化を高速に処理可能な並列処理方法に関する。
【背景技術】
【０００２】
実対称行列の固有地問題を解く上で対角化を行い、三重対角行列の固有値問題に変形して解く方法がある。この方法では、実対称行列を三重対角化する部分の計算コストが大きい。この三重対角化に関して、共有メモリ型の並列計算機では、ベクトル計算機のような強力なメモリアクセス能力がないため、メモリアクセスに対して演算量を増やすアルゴリズム上の工夫がなされている。【０００３】
特許文献１には、ブロック化された三重対角化の処理の基本が記述されている。
しかし、計算のコストの大部分が行列ベクトル積であるため、メモリアクセスのスピードの影響を大きく受ける。全体の性能を改善するために、この部分の改良が必要である。【０００４】
演算量を増やす工夫を行ったブロック化した方法では、行列ベクトル積を連続アクセスを使って計算すること、および、その計算を各ＣＰＵで均等に並列処理することがポイントになる。このために、更新時に必要となる行列全体を更新するとき、下三角部分を更新し、上三角部分をコピーしていた。更新時の負荷をバランスさせるために、行列を列方向に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:An のn+i 行ベクトルよりHouseholder ベクトルuを作る。
step2: vi = An +iu , wi = vi - u(utv)/2, を計算する。
step3: Ui =(Ui-1,ui ) , Wi =(Wi-1,wi )と更新する。
step4: if(i < blks)then 
Anのn+i+1 列目を更新する。【００１３】
An(*,n+i+1) = An (*,n+i+1)- Ui Wi (n+i+1,*)t- Wi Ui (n+i+1,*)t 
endif
enddo
step５：An+blks = An - Ublks Wblkst- Wblks Ublkst
次に、ブロック化されたHouseholder 変換による三重対角化について説明する。ブロック化されたHouseholder変換は、以下の式の通りに行う。
ｖ＝（v1,v2,...,vn )     
｜ｖ｜2 ＝ｖ・ｖ＝ｈ2 、
Uv =(h,0,…,0)t とすると、Uv =  v  −  （v1−h ，v2，．．．，vn）の関係がある。
U ＝（I −uut / ｜u ｜2 ）＝（I −αuut ）, ここでu=(v1−h ，v2，．．．，vn ) 
以下αを省略して計算する。
An+1 ＝Ut An U ＝（I-uut ）A （I-uut）
＝An - uut An  -  An uut＋  uut An uut 
＝An - uvt  -vut ＋uvt uut
＝An - uwt -uut ut v/2 - wut- uut ut v/2 +uut ut v
＝An - u wt  -wut       ・・・（＊）
ここで、  w = v - u(ut v)/2,  v = An u 
これを繰りし利用すれば、
An+k   ＝An   −Uk Wkt   −WkUkt     ・・・（＊＊）
ｋステップ目の計算は（＊）および（＊＊）からVn を次のように計算できる。
vk = An uk -Uk-1 Wk-1tuk - Wk-1Uk-1t uk 
wk =vk - uk ukt vk/2 
Uk =(Uk-1,uk )  ,  Wk =(Wk-1,wk) 
An+k = An - Uk Wkt - WkUkt 
ｂ）Householder変換を構成する情報の格納
固有ベクトルを計算するときに３重対角化でつかったHouseholder 変換が必要になる。このため、Ｕn およびαをhouseholder 変換を構成するベクトルの位置に格納する。αは対応する対角要素の位置に格納する。
ｃ）Ｕi を効率的に求める方法
ブロック部分の三重対角化を行うためのHouseholder 変換のために、次のベクトルを更新する必要がある。これらの計算をできるだけローカルに行うために、ブロック幅部分を作業域にコピーして３重対角化を行い、もとの領域に格納する。次の列ベクトルの更新をそのたび行わず、行列積の形にして演算密度を上げて計算を行う。このため、ブロック部分の三重対角化を再帰的なプログラミングで行う。【００１４】
すなわち、図３の左の図における行列Ａのブロック領域を、図３右の作業領域Ｕにコピーし、当該ブロック領域を前半更新部分と後半更新部分とに分けて、再帰的に処理する。図４Ａ〜図４Ｃは、再帰的な処理の様子を示している。今、図４Ａの斜線部分を更新したとすると、次に、図４Ｂの斜線部分を更新し、更に、図４Ｃの斜線部分を更新する。このように、前半部分と後半部分とに分けて更新する仕方を入れ子状に適用する。すなわち、再帰的プログラムが深さ２まで呼ばれると、最初の前半処理で図４Ａの斜線部分をＢとして更新し、次に、図４Ｂの斜線部分、最後に図４Ｃの斜線部分が更新される。更新時の並列化は更新部分ブロック行列を行ベクトル方向に均等に分割して、各部分の更新をプロセッサで均等に分割した領域に対して並列して行う。【００１５】
図５〜図７は、本発明の実施形態に従った並列処理方法を説明する図である。
１）行列ベクトル積を下三角部分のみを使って１つのアクセス（ロード）で 2倍の演算を行ない、各ＣＰＵで均等に行なう方法
図５に示されるように、下三角行列部分のみを利用して行列ベクトル積を計算する。行列ベクトル積で使う行列の大きさをＮｎ×Ｎｎとする。４並列で行うとき、下三角部分を４つの部分に分割する位置は、以下のようになる。
im=Nn×(1-sqrt(1-m/#p))
ｍ＝１、２、３で計算したものがちょうど境界位置になる。なお、＃ｐは、分割する総数、つまり、並列数である。【００１６】
Ａnの対称性つまり、Ａn ＝ Ａnt を利用して、各プロセッサで自分の受け持つ領域にある情報から、行列積を計算する。
各部分は対角ブロック行列部分の下三角行列部分Ｄm と、その下にある長方形の行列Ｒmとからなる。【００１７】
下三角行列を#Pに分割した、m 番目の部分を転置したものは、ちょうど上三角部分を行方向に分割したものになっている。この関係から、以下の計算を行なうことができる。
ｖxm1 (im＋１：Ｎn )=Ｒm ×ｕ1(ｉm-1＋１：ｉm)           （式１）
ｖxm2 (im-1＋１：ｉm )=Ｒmt ×ｕ1(ｉm＋１：Ｎn)           （式２）
ｖxm3 (im-1＋１：ｉm )=Ｄm ×ｕ1(ｉm-1＋１：ｉm)           （式3 ）
ｖxm4 (im-1＋１：ｉm )=ｎｏｔｄｉａｇ（Ｄm ）t×ｕ1(ｉm-1＋２：ｉm ) （式4 ）
（ｎｏｔｄｉａｇ（Ｘ）は、下三角行列Ｘの対角部分を除いた下三角行列を表す）【００１８】
【数１】
2008136045.tif 000003
この計算で、行列の要素を連続にアクセスしながら、１つの要素のアクセスに対して、ｖ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)のアクセスに比べ、演算はほとんど無視できるほど速いので、別々に行列ベクトル積を行なう場合に比べて約２倍の性能になる。【００２１】
（＊）分割の境界を計算する式、im = Nn ×(1-qsrt(1-m/#p))は次のように求めることができる。
各CPU での演算数が等しい。つまり、要素数に比例する行列の対応部分の面積（あるいは要素数そのもの）が等しいと考える。
Ｎn2-( Ｎn - ｉm )2＝Ｎn2×m/#p，ｒ＝ｉm ／Ｎn と置いて計算する。【００２２】
ｒ2 −２ｒ＝m/#p    ，  この根の0<= ｒ<=1 をとる。
(ｒ−１)2＝m/#p    
ｒ＝１−sqrt(1-m/#p)  ( 根の範囲より小さい方を採用する。）
２）更新部分の並列化および対角ブロック部分の計算方法
step5の更新部分は対称性を利用して下半分を計算する。更新の演算量が各CPU(スレッド) で、均等化するように、分割する。分割方法は行列ベクトル積で分割したのと同じ方法で行なう。【００２３】
行列ベクトル積で使う行列の大きさをNn ×Nn とする。4 並列で行なうとき、下三角部分を 4つの部分に分割する位置は以下のようになる。
i m = Nn ×(1-qsrt(1-m/#p))      m=1,2,3 で計算したものがちょうど境界の位置になる。#pは、分割する総数、つま並列数である。【００２４】
このように、ＣＰＵごとに分割された領域は、対角ブロック行列部分とその下の長方行列部分に分かれる。図６は、配列An+iをＣＰＵごとに分割した様子を示している。長方行列部分は、普通の行列積で
An+k ＝An −Uk Wkt −WkUkt から、対応部分を計算する。
対角ブロック行列部分の更新は、図７のように行なう。すなわち、対角ブロックの下三角部分について、対角線の２等分点を決める。この等分点から、行方向、及び、列方向に下三角行列を分割する。矩形部分は、行列積で更新する。２つの大きさが半分の下三角行列ができるので、これらに関して、同じく対角線上の等分点を決め、同様に列方向、行方向に分割する。そして、対象となる下三角行列の大きさが十分小さくなったら、上三角行列も一緒に行列積で計算する。これらは、簡単な再帰プログラムで作ることが出来る。【００２５】
以上の並列処理方法により、高性能かつスケーラビリティのある実対称行列の三重対角化を実現できる。以上の方法は、従来の方法に比べて２倍超の性能であることがわかった。【００２６】
図８〜図１１は、本発明の実施形態の並列処理方法をより詳しく説明する図である。
ブロック化された三重対角化演算においては、以下のような処理を行う。
１）三重対角化でブロックごとに、三重対角化を行う。上のブロックを三重対角化し、その三重対角化で生成された行列U およびW を使って、ブロック幅分縮小された正方行列( この下三角行列は、図８において、太い点線で表示) を更新する。ただし、この正方行列は対称行列なので、下三角行列部分( 太い点線で囲んだ三角形部分) のみ更新して上三角行列部分のデータとして利用する。
２）ブロックごとの三重対角化では、ブロック内の太線で示した正方行列(An+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)=L2t*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である。
ステップＳ７０において、バリア同期を取り、ステップＳ７１において、【００４４】
【数２】
2008136045.tif 000004
を計算し、このサブルーチンを抜ける。
図１７は、サブルーチン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*Ut)を行う。サブルーチンltgemmtrtを呼び出し、対角ブロック行列の下三角行列部分を更新する。引数として、W(is:ie, 1:blk), U(is:ie, 1:blk), A(is:ie, is:ie)及び対角ブロック行列の大きさlen=ie-is+1を渡す。【００４８】
ステップＳ９９において、下三角行列部の更新(A=A-U*Wt)を行う。サブルーチン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を渡す。そして、ステップＳ１１３の後、このサブルーチンを抜ける。
【図１】
2008136045.tif 000005
【図２】
2008136045.tif 000006
【図３】
2008136045.tif 000007
【図４Ａ】
2008136045.tif 000008
【図４Ｂ】
2008136045.tif 000009
【図４Ｃ】
2008136045.tif 000010
【図５】
2008136045.tif 000011
【図６】
2008136045.tif 000012
【図７】
2008136045.tif 000013
【図８】
2008136045.tif 000014
【図９】
2008136045.tif 000015
【図１０】
2008136045.tif 000016
【図１１】
2008136045.tif 000017
【図１２】
2008136045.tif 000018
【図１３】
2008136045.tif 000019
【図１４】
2008136045.tif 000020
【図１５】
2008136045.tif 000021
【図１６】
2008136045.tif 000022
【図１７】
2008136045.tif 000023
【図１８】
2008136045.tif 000024
【図１９】
2008136045.tif 000025
【図２０】
2008136045.tif 000026
【国際調査報告】
2008136045.tif 000027
2008136045.tif 000028
2008136045.tif 000029