On 18-06-2013, at 13:23, peter dalgaard <pda...@gmail.com> wrote:
> On some further digging, on the Lion machine, the problem seems absent from
> builds against the veclib framework. I strongly suspect that the issue is the
> same as PR#14964, which also affects Hermitian matrices of dimension 33x33
> and above.
>
With the CRAN distribution of R-3.0.1 on OS X 10.8.4 I encountered the
problem too.
I've done a small experiment to see if I could find a possible cause.
This version of R uses the libRblas and libRlapack.
For complex matrices I'm assuming that eigen uses the Lapack routine zheev and
that R first calls zheev to determine the optimal workspace.
I've made an alternative eigen that also uses Lapack's zheev but calls it with
minimal workspace to force zheev to use the unblocked algorithm. To determine
if the blocked algorithm that Lapack uses could be the cause of the problem.
This will work when eigen is called first so that the Lapack library is
dyn.loaded.
Code:
# quick and dirty
zeigen <- function(A,symmetric,only.values=FALSE) {
# SUBROUTINE ZHEEV( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, RWORK,
# * INFO )
# .. Scalar Arguments ..
# * CHARACTER JOBZ, UPLO
# * INTEGER INFO, LDA, LWORK, N
# * ..
# * .. Array Arguments ..
# * DOUBLE PRECISION RWORK( * ), W( * )
# * COMPLEX*16 A( LDA, * ), WORK( * )
n <- dim(A)[1]
# use minimal workspace
lwork <- as.integer(2*n-1)
# warning: "N" and "L" work on OSX but you may have to write an interface
routine since some systems
# give Fortran runtime errors when passing character arguments to Fortran
routines.
z <- .Fortran("zheev", "N", "L", n, A=A, n, w=numeric(n),
work=complex(lwork), lwork,
numeric(3*n-2), info=integer(1L))
list(values=z$w, vectors=z$A)
}
N <- 100
jj <- matrix(0,N,N)
A <- exp(-0.1*(row(jj)-col(jj))^2)
min(eigen(A,only.values=TRUE,symmetric=TRUE)$values)
## [1] 2.521153e-10
# Coercing A to a complex matrix should make no difference, but makes
# eigen() return the wrong answer:
min(eigen(A+0i,only.values=TRUE,symmetric=TRUE)$values)
## [1] -0.359347
min(zeigen(A+0i,only.values=TRUE,symmetric=TRUE)$values)
##[1] 2.521154e-10
So it seems that the blocked algorithm is the cause of the error and that using
the (possibly slow) unblocked algorithm gives the correct result.
Berend
______________________________________________
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
