From ripley at stats.ox.ac.uk  Fri Nov  5 23:10:19 2004
From: ripley at stats.ox.ac.uk (Prof Brian Ripley)
Date: Fri Nov  5 23:10:24 2004
Subject: [Rd] Re: [R] fast partial spectral decompositions.
In-Reply-To: <2C0E69E5-2F6A-11D9-AC45-000A95DA144C@mcauliffe.com>
Message-ID: <Pine.LNX.4.44.0411052110270.25202-100000@gannet.stats>

[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


[Rd] Re: [R] fast partial spectral decompositions.
Prof Brian Ripley