matrix - BLAS 矩阵乘矩阵转置乘法

Question

我必须以 A'A 的形式计算一些产品或更一般的 A'DA ，其中 A是将军mxn矩阵和 D是对角线 mxm矩阵。他们都是满级；即 rank(A)=min(m,n) .

我知道你可以节省大量时间是这样的对称产品:鉴于A'A是对称的，您只需要计算乘积矩阵的下-或上-对角线部分。这增加了 n(n+1)/2要计算的条目，大约是典型 n^2 的一半对于大型矩阵。

这是我想利用的一个很好的节省，而且我知道我可以在 for 中实现矩阵-矩阵乘法。环形 。但是，到目前为止，我一直在使用 BLAS，它比任何 for 都要快得多。我可以自己编写的循环实现，因为它优化了缓存和内存管理。

有没有办法有效地计算A'A甚至 A'DA使用BLAS？
谢谢!

Answer 1

您正在寻找 dsyrk BLAS 的子程序。

如文档中所述:

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

DSYRK performs one of the symmetric rank k operations

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

or

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

where alpha and beta are scalars, C is an n by n symmetric matrix and A is an n by k matrix in the first case and a k by n matrix in the second case.

在 A'A的情况下存储上三角是:

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

对于 A'DA BLAS 中没有直接的等价物。但是您可以使用 dsyr在 for 循环中。

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

DSYR performs the symmetric rank 1 operation

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

where alpha is a real scalar, x is an n element vector and A is an n by n symmetric matrix.

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