LAPACKer's,
Over the years I've repeatedly run into a symmetric matrix modification
operation that's not captured in the BLAS. Following Beresford's
symmetric letter convention, the operation is
A = A + C^T W C
where A and W are real symmetric (or complex Hermitian) and C is
rectangular. This occurs as a rank 2 modification in symmetric
indefinite factorizations. You've heard me before on that. Now I've
encountered it in a very general form in computational chemistry. Here C
and W are full rank matrices, of considerably large rank.
Obviously this could be done with one intermediate step, as A + (C^T M)
C or A + C^T (M C). Storing the intermediate product creates a storage
penalty that's not particularly important at the moment. The major
problem is that the only place where BLAS or the PBLAS recognize that a
matrix-matrix product X*Y is symmetric is when Y = X^T. (The symmetric
matrix argument to PDSYMM is X or Y, but the product is general.) This
means that we have only PDGEMM available to do this, which requires
twice as much work as we really need to do.
Has anyone else ever faced this problem? Does anyone know of an
implementation of either
* the elegant solution (modifying PDSYRK to accept the extra
symmetric matrix argument and to do the intermediate products on
the fly), or
* the inelegant solution, computing the intermediate product as a
separate matrix and then modifying PDGEMM to avoid the work in one
or the other triangle?
Thanks!
John
|