使用 gsl、Blas 和 Lapack 计算 (Aᵀ×A)*(Bᵀ×B) 矩阵答案

【问题标题】：Compute (Aᵀ×A)*(Bᵀ×B) matrix using gsl, Blas and Lapack使用 gsl、Blas 和 Lapack 计算 (Aᵀ×A)*(Bᵀ×B) 矩阵
【发布时间】：2019-02-14 05:00:51
【问题描述】：

让我有 2 个对称矩阵：

A = {{1,2}, {2,3}}
B = {{2,3},{3,4}}

我可以使用 gsl、Blas 和 Lapack 计算矩阵 (A^T×A)*(B^T×B) 吗？

我正在使用

gsl_blas_dsyrk(CblasUpper, CblasTrans, 1.0, A, 0.0, ATA);
gsl_blas_dsyrk(CblasUpper, CblasTrans, 1.0, B, 0.0, BTB);
gsl_blas_dsymm(CblasLeft, CblasUpper, 1.0, ATA, BTB, 0.0, ATABTB); // It doesn't work

(Aᵀ·A) = ATA = {{5, 8}, {0, 13}} -- ok, gsl_blas_dsyrk returns symmetric matrix as upper triangular matrix.
(Bᵀ·B) = BTB = {{13, 8}, {0, 25}} -- ok.
(Aᵀ·A)·(Bᵀ·B) = ATABTB = {{65, 290}, {104, 469}} -- it's wrong.

【问题讨论】：

标签： matrix matrix-multiplication lapack blas gsl

【解决方案1】：

对称BTB，问题就解决了。

如您所见，对称矩阵的上三角部分由dsyrk() 计算。然后应用dsymm()。根据dsymm()的定义，由于使用了CblasLeft标志位，所以进行如下操作：

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

其中 alpha 和 beta 是标量，A 是对称矩阵，B 和 C 是 m × n 矩阵。

确实，B 矩阵是一般矩阵，不一定是对称矩阵。结果，ATA 乘以 BTB 的上三角部分，因为 BTB 的下三角部分不计算。

对称BTB，问题就会解决。为此，for循环是一个简单的解决方案，请参阅Convert symmetric matrix between packed and full storage?

【讨论】：