BLAS Matrix durch transponierte Matrix Multiply

Question

I einige Produkte in der Form zu berechnen haben A'A oder allgemeinere A'DA, wo A eine allgemeine mxn Matrix ist und D ist eine diagonale Matrix mxm. Beide sind voller Rang; d.h. rank(A)=min(m,n).

Ich weiß, dass Sie eine erhebliche Zeit sparen können, ist solche symmetrische Produkte: vorausgesetzt, dass A'A symmetrisch ist, müssen Sie nur den unteren - oder oberen - diagonalen Teil der Produktmatrix berechnen. Dies addiert sich zu n(n+1)/2 zu berechnenden Einträgen, was in etwa der Hälfte des typischen n^2 für große Matrizen entspricht.

Dies ist eine große Einsparung, die ich ausnutzen möchte, und ich weiß, dass ich die Matrix-Matrix-Multiplikation innerhalb einer for Schleife implementieren kann. Bisher habe ich jedoch BLAS verwendet, was viel schneller ist als jede andere for Schleifenimplementierung, die ich selbst schreiben könnte, da sie die Cache- und Speicherverwaltung optimiert.

Gibt es eine Möglichkeit, A'A oder sogar A'DA mithilfe von BLAS effizient zu berechnen? Danke!

Answer 1

Sie suchen nach Unterprogramm dsyrk von BLAS.

Wie in der Literatur beschrieben:

SUBROUTINE dsyrk (Uplo, TRANS, N, K, ALPHA, A, LDA, BETA, C, LDC)

DSYRK führt eine der symmetrischen rank k Operationen

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

oder

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

wo Alpha und Beta sind Skalare, C ist eine n von n symmetrische Matrix und A ist eine n von k-Matrix im ersten Fall und eine k-n-Matrix im zweiten Fall.

Im Fall von A'A oberes Dreieck Speicherung:

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

Für die A'DA es keine direkte Entsprechung in BLAS ist. Sie können jedoch dsyr in einer for-Schleife verwenden.

SUBROUTINE dsyr (Uplo, N, ALPHA, X, INCX, A, LDA)

DSYR führt den symmetrischen Rang 1 Betrieb

A := alpha*x*x**T + A,

wobei alpha ein echter Skalar , x ist ein n-Element-Vektor und A ist eine n mal n-symmetrische Matrix.

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

Answer 2

SYRK ist für A'A in Ordnung. Für A'DA können Sie SYMM für eine Seite verwenden, zB V = A'D und dann verwenden Sie intel MKL's GEMMT für W = V A.GEMMT ist wie GEMM, außer dass es die Tatsache ausnutzt, dass die Ergebnismatrix symmetrisch ist und somit nur etwa die Hälfte der Arbeit erledigen muss.