[R] Understanding dsyrk_ in C code

[Douglas Bates]

On Tue, Jan 6, 2009 at 8:54 PM, Nathan S. Watson-Haigh
<nathan.watson-haigh at csiro.au> wrote:
> -----BEGIN PGP SIGNED MESSAGE-----
> Hash: SHA1
>
> I'm trying to understand some C code in an R package I'm using. I'm address this question here as
> it's matrix algebra...and I'm no pro at that!
>
> the C command reads:
>
> double alpha = 1.0, beta = 0.0;
> dsyrk_("L", "N", nGenes, nGenes, & alpha, mat1, nGenes,
>         & beta, mat2, nGenes);
>
> - From google, I've found out that dsyrk is for performing one of the symmetric rank k operations -
> whatever that means!? From here:
> http://linux.die.net/man/l/dsyrk
>
> I've found that the calculation being performed is:
> alpha*A*A' + beta*C
>
> However, since alpha is 1 and beta is 0, this reduces to:
> => 1*A*A' + 0*C
> => A*A'
>
> Which is simply the cross product....am I correct?

Not quite.  The crossprod function in R computes A'A.  This particular
operation is what is computed by the tcrossprod function in R.

Another difference is that this call modifies only the diagonal and
lower triangle of mat2.  The tcrossprod function calls the same
underlying code then copies the lower triangle to the upper triangle
so as to produce a symmetric matrix.

When you are calling Lapack or BLAS routines from C it is helpful to
look at the declarations in include/R_ext/BLAS.h and
include/R_ext/Lapack.h

/* DSYRK - perform one of the symmetric rank k operations */
/* C := alpha*A*A' + beta*C or C := alpha*A'*A + beta*C */
BLAS_extern void
F77_NAME(dsyrk)(const char *uplo, const char *trans,
		const int *n, const int *k,
		const double *alpha, const double *a, const int *lda,
		const double *beta, double *c, const int *ldc);

This helps to explain the mysterious "L" and "N" as the first two
arguments in the call.

By the way, the use of the literal name dsyrk_ in that code is risky.
It is preferable to use the macro F77_SUB(dsyrk) which is defined in
the include file include/R_ext/RS.h

If you are wondering what the difference between F77_SUB and F77_NAME
might be, there isn't a difference in R.  Long, long ago there was a
difference between the two macros in S and S-PLUS related to a
peculiar Fortran compiler, if I recall correctly.  Most of the time
these macros simply append an underscore to the name but there may
still be cases where something else occurs and it is best to allow for
those cases.