BLAS矩阵乘矩阵转置乘法答案

【问题标题】：BLAS matrix by matrix transpose multiplyBLAS矩阵乘矩阵转置乘法
【发布时间】：2018-04-11 08:01:11
【问题描述】：

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

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

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

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

【问题讨论】：

标签： matrix linear-algebra blas

【解决方案1】：

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

如文档中所述：

子程序 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 矩阵，在第二种情况下是 k × n 矩阵。

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

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

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

子程序 dsyr(UPLO,N,ALPHA,X,INCX,A,LDA)

DSYR 执行对称秩 1 操作

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

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

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

【讨论】：

【解决方案2】：

SYRK 适合 A'A。对于 A'DA，您可以将 SYMM 用于其一侧，例如 V = A'D，然后将 intel MKL 的 GEMMT 用于 W = V A。GEMMT 与 GEMM 类似，只是它利用了结果矩阵是对称的这一事实，并且因此只需要完成大约一半的工作。

【讨论】：

【解决方案3】：

@ztik 建议的 BLAS 中的 dsyrk 例程是 A'A 的例程。对于A'DA，一种可能性是使用可以执行对称秩2k 操作的dsyr2k 例程：

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

设置alpha = 0.5, beta = 0.0，并设置B = DA。请注意，这种方式假定您的对角矩阵 D 是真实的。

【讨论】：