2017-10-30 13 views
2

IはA一般mxn行列であり、Dは対角行列であるmxm形態A'A又はA'DAより一般的にいくつかの製品を計算する必要が転置。どちらもフルランクです。すなわちrank(A)=min(m,n)BLAS行列が乗算

私はあなたが相当な時間を節約することができることを知っています。そのような対称的な製品です:A'Aが対称であるとすれば、製品マトリックスの下側 - または上側 - 対角線部分を計算すればよい。すなわち、概ね大きな行列のための典型的なn^2の半分である、計算するn(n+1)/2エントリに追加します。

これは私が活用したい素晴らしいセービングである、と私は私がforループ内で行列 - 行列乗算を実装することができます知っています。しかし、これまでのところ、私はそれをキャッシュとメモリ管理を最適化するので、私は、自分で書くことができます任意のforループの実装よりもはるかに高速であるBLASを、使用しています。

効率的にBLASを使用してA'AかさえA'DAを計算する方法はありますか? ありがとう!

答えて

2

あなたはBLASのdsyrkサブルーチンのために見ています。

ドキュメントに記載されているように:

dsyrkの対称ランクのいずれかを行う

SUBROUTINEのdsyrkの(UPLO、TRANS、N、K、ALPHA、A、LDA、BETA、C、LDC) K動作

C := alpha*A*A**T + beta*C

又は

C := alpha*A**T*A + beta*Cは、αおよびβはスカラーである

は、Cは、n行n列の対称マトリックスであり、Aは、最初のケースにおけるk行列によってNと第2ケースにおけるn行列によってKです。上三角を記憶A'Aの場合

である:A'DAについて

CALL dsyrk('U' , 'T' , N , M , 1.0 , A , M , 0.0 , C , N) 

BLASには直接等価物は存在しません。しかし、forループでdsyrを使用することができます。

SUBROUTINEのdsyr(UPLO、N、ALPHA、X、INCX、A、LDA)

DSYRは、対称ランク1操作

A := alpha*x*x**T + Aを行うアルファは実数スカラー

xはn要素ベクトル、Aはn×n対称行列です。

do i = 1, M 
    call dsyr('U',N,D(i,i),A(1,i),M,C,N) 
end do 
+0

素晴らしいですね。ありがとう! – enanone

-1

SYRKはA'Aの罰金です。 A'DAの場合、SYMMを使用することができます。たとえば、V = A'Dのようにして、W = V Aの場合はintel MKLのGEMMTを使用します。GEMMTはGEMMと似ていますが、結果行列が対称であることを利用する点を除いて、ほぼ半分の作業しか必要としません。