Efficient way of computing X'*A*X where A is symmetric

[eev234]

Hi All.

I was would like to ask what is the best way of evaluating the following product of matrices:

Let X be m-by-n, A be m-by-m and symmetric and I want to compute
C = X' * A * X
where X' is the transpose of X. The result is a symmetric n-by-n matrix. The only option I can find is to call dgemm twice. Any ideas?

[cottrell]

So far as I know, there is no built-in BLAS function that will do this.

You could take a look at the implementation in libgretl, at
http://gretl.cvs.sourceforge.net/viewvc/gretl/gretl/lib/src/gretl_matrix.c?revision=1.495&view=markup
Look for gretl_matrix_qform (for quadratic form).

The "gretl_matrix_set" and "gretl_matrix_get" things are macros. defined in
http://gretl.cvs.sourceforge.net/viewvc/gretl/gretl/lib/src/gretl_matrix.h?revision=1.237&view=markup

The nice thing about the "gretl_matrix" type from the pov of this forum
is that the actual data storage is just a one-dimensional array of doubles,
in column-major order, and so is fully compatible with lapack.

[eev234]

Dear cottrell,

Thank you for your reply. That was very helpful. I looked into it and I think it's what I am looking for.