Compute (Aᵀ×A)*(Bᵀ×B) matrix using gsl, Blas and Lapack

Question

Let I have 2 symmetric matrices:

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

Can I compute the matrix (A^T×A)*(B^T×B) using gsl, Blas and Lapack?

I'm using

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

It returns:

(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.

Answer 1

Symmetrize BTB and the problem will be solved.

As you noticed, the upper triangular parts of symmetric matrices are computed by dsyrk(). Then dsymm() is applied. According to the definition of dsymm(), the following operation is performed since the flag CblasLeft is used:

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

where alpha and beta are scalars, A is a symmetric matrix and B and C are m by n matrices.

Indeed, the B matrix is a general matrix, not necessarly a symmetric one. As a result, ATA is multiplied by the upper triangular part of BTB, since the lower triangular part of BTB is not computed.

Symmetrize BTB and the problem will be solved. To do so, for loops is a straightforward solution , see Convert symmetric matrix between packed and full storage?