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

您正在寻找dsyrkBLAS的子程序。

如文档中所述：

SUBROUTINE dsyrk（UPLO，TRANS，N，K，ALPHA，A，LDA，BETA，C，LDC）

DSYRK执行对称秩k运算之一

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

要么

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

其中alpha和beta是标量，C在第一种情况下是n×n对称矩阵，A在第一种情况下是n×k矩阵，第二种情况是ak×n矩阵。

在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执行对称等级1操作

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

其中alpha是实数标量，x是n个元素向量，A是n x n对称矩阵。

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