R help - Nov 2004 - fast partial spectral decompositions.

hello,

i want to compute the top k eigenvalues+eigenvectors of a (large)
real symmetric matrix. since it doesn't look like any top-level R
function does this, i'll call LAPACK from a C shlib and then
use .Call. the only LAPACK function i see to do this in
R_ext/Lapack.h is dsyevx. however, i know that in LAPACK dsyevr
can also return a partial eigendecomposition. why is dsyevr not
exported in R_ext/Lapack.h? my superficial understanding is that
dsyevr is "better" (faster? stabler?) for both complete and
partial eigenproblems than dsyevd/dsyevx, but only the complete
eigenproblem interface to dsyevr appears to be exported in
Lapack.h (as dsyev).

corrections to misunderstandings in the above are welcome. advice
on whether using dsyevr rather than dsyevx is (very) important
for partial decompositions is also gratefully accepted.

please include jon at mcauliffe.com in the reply.

jon.

[This really is a programming question which the posting guide says should
have been sent to R-devel, so I have diverted it there.]

dseyv is not an interface to dsyevr, but a separate routine.

R does use dsyevr these days, but before R required IEC60559 arithmetic, 
it also provided the choice of dsyev, as the latter does not require 
IEC60559.

I am not sure what R_exts/Lapack.h is intended to be.  If you include
$(LAPACK_LIBS) $(BLAS_LIBS) in the building of your package, you will get
access to a full double-precision LAPACK library including dsyevr.
However, that header file is _not_ part of the R API and is _not_ a list
of exports.

There is another complication.  R imports LAPACK subroutines from either
an external LAPACK library or from an LAPACK library it creates.  Because
some external LAPACK libraries have a broken dsyev(r) but might be in use
as a BLAS library, under some circumstances R uses a renamed dsyev(r) as
rsyev(r).  So it is potentially dangerous to make use of dsyev(r) (and, 
let me say again, they are not part of the R API).  (AFAIR we did this to 
avoid getting incorrect results on some version of libsunperf.)

I suggest you use your own LAPACK routines of known provenance.

On Fri, 5 Nov 2004, Jon McAuliffe [malfunctioning@shift.key] wrote:
> i want to compute the top k eigenvalues+eigenvectors of a (large)
> real symmetric matrix. since it doesn't look like any top-level R
> function does this, i'll call LAPACK from a C shlib and then
> use .Call. the only LAPACK function i see to do this in
> R_ext/Lapack.h is dsyevx. however, i know that in LAPACK dsyevr
> can also return a partial eigendecomposition. why is dsyevr not
> exported in R_ext/Lapack.h? my superficial understanding is that
> dsyevr is "better" (faster? stabler?) for both complete and
> partial eigenproblems than dsyevd/dsyevx, but only the complete
> eigenproblem interface to dsyevr appears to be exported in
> Lapack.h (as dsyev).
> 
> corrections to misunderstandings in the above are welcome. advice
> on whether using dsyevr rather than dsyevx is (very) important
> for partial decompositions is also gratefully accepted.
I am baffled by the problem with your shift key: it does work some of the
time but your text is very hard to read when it does not work.

-- 
Brian D. Ripley,                  ripley@stats.ox.ac.uk
Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
University of Oxford,             Tel:  +44 1865 272861 (self)
1 South Parks Road,                     +44 1865 272866 (PA)
Oxford OX1 3TG, UK                Fax:  +44 1865 272595