I'm trying to solve the linear system AXA=B with LAPACK. Matrices A and B are complex and symmetric 1000x1000. The matrix B is 70% filled with zeros and block diagonal but A is dense.

PROCEDURE 1

To solve the problem I rewrite the problem as AY=B, solves for Y with zgetrs and zgetrf.

After all, solves the system Y=XA as (XA)^T = Y^T = A^T*X^T => AX^T = Y^T.

PROCEDURE 2

To confirm that result, I use two dgemm calls to regain AXA that I expected it to be the matrix B.

The problem: The result of PROCEDURE2 don't agree with B! the Frobenius norm (zlange) of AXA after PROCEDURE2 is one order of magnitute grater than the original B matrix and is not at all symmetric!

I really appreciate suggestions! as don't understand what is going wrong with this very simple (5 LAPACK calls) procedure.

Thanks in advanced!