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=markupLook 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=markupThe 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.