Re: [Rd] eigen(symmetric=TRUE) for complex matrices

2013-06-19 Thread peter dalgaard
On Jun 18, 2013, at 21:49 , peter dalgaard wrote: > > On Jun 18, 2013, at 21:23 , Berend Hasselman wrote: > >> >> 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. > > Thanks Berend, > > The f

Re: [Rd] eigen(symmetric=TRUE) for complex matrices

2013-06-19 Thread Berend Hasselman
On 19-06-2013, at 10:24, peter dalgaard wrote: > > On Jun 18, 2013, at 21:49 , peter dalgaard wrote: > >> >> On Jun 18, 2013, at 21:23 , Berend Hasselman wrote: >> >>> >>> So it seems that the blocked algorithm is the cause of the error and that >>> using the (possibly slow) unblocked algo

Re: [Rd] eigen(symmetric=TRUE) for complex matrices

2013-06-19 Thread Ravi Varadhan
I am also able to reproduce this problem on Windows: library(eiginv) n <- 33 # problem does not arise for n <= 32 evals <- sort(rnorm(n)) system.time(A <- eiginv(evals, symmetric=TRUE)) all.equal(evals, sort(eigen(A+0i,,TRUE)$val)) cbind(evals, sort(eigen(A+0i,,TRUE)$val)) Best, Ravi >

Re: [Rd] eigen(symmetric=TRUE) for complex matrices

2013-06-19 Thread peter dalgaard
On Jun 19, 2013, at 12:43 , Berend Hasselman wrote: > > On 19-06-2013, at 10:24, peter dalgaard wrote: > >> >> On Jun 18, 2013, at 21:49 , peter dalgaard wrote: >> >>> >>> On Jun 18, 2013, at 21:23 , Berend Hasselman wrote: >>> So it seems that the blocked algorithm is the cause

Re: [Rd] eigen(symmetric=TRUE) for complex matrices

2013-06-19 Thread Roger Bivand
peter dalgaard gmail.com> writes: > > ... > > There could be compiler issues. It could also be that the internal LAPACK uses a newer version than the > external libs and a bug was introduced in between them. Running: N <- 100 jj <- matrix(0,N,N) A <- exp(-0.1*(row(jj)-col(jj))^2) min(eigen(A

Re: [Rd] eigen(symmetric=TRUE) for complex matrices

2013-06-19 Thread Berend Hasselman
On 19-06-2013, at 14:17, peter dalgaard wrote: > > > > Thanks. I think I have it nailed down now. The culprit was indeed in our > reference BLAS (I had only checked the LAPACK code), cmplxblas.f to be > specific. Revision 53001 had a number of IF statements being commented out, > but two o

Re: [Rd] eigen(symmetric=TRUE) for complex matrices

2013-06-19 Thread peter dalgaard
On Jun 19, 2013, at 15:58 , Berend Hasselman wrote: > > On 19-06-2013, at 14:17, peter dalgaard wrote: > >> >> >> >> Thanks. I think I have it nailed down now. The culprit was indeed in our >> reference BLAS (I had only checked the LAPACK code), cmplxblas.f to be >> specific. Revision 530

[Rd] Calling an array in a struct in C to R

2013-06-19 Thread Tee-Jay-Ardie
Hi there, Although I'm a quite experienced R user and know my ways in C, I stumbled upon a problem I don't know how to solve. Therefore, I hope someone can provide me with the information or pointers I need in order to understand the way in which the communication between R and C occurs. I have th

Re: [Rd] Calling an array in a struct in C to R

2013-06-19 Thread William Dunlap
> should I rewrite my C code using Call in R? Yes, rewrite it to use the .Call() interface. Bill Dunlap Spotfire, TIBCO Software wdunlap tibco.com > -Original Message- > From: r-devel-boun...@r-project.org [mailto:r-devel-boun...@r-project.org] On > Behalf > Of Tee-Jay-Ardie > Sent: We

[Rd] str and s3 classes

2013-06-19 Thread Hadley Wickham
Hi all, Because str uses the generic version of length and names, it's currently very easy to create objects that break str: a <- structure(list(1:5), class = "a") length.a <- function(x) 2L str(a) b <- structure(list(1:5), class = "b") names.b <- function(x) character() str(b) These are obv

Re: [Rd] str and s3 classes

2013-06-19 Thread Duncan Murdoch
On 13-06-19 4:44 PM, Hadley Wickham wrote: Hi all, Because str uses the generic version of length and names, it's currently very easy to create objects that break str: a <- structure(list(1:5), class = "a") length.a <- function(x) 2L str(a) b <- structure(list(1:5), class = "b") names.b <- fu