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

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 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 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 Ravi Varadhan
ct.org Subject: Re: [Rd] eigen(symmetric=TRUE) for complex matrices 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 33x

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 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-18 Thread peter dalgaard
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 finger seems to be pointing to the internal libRlapack/libRblas. The

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

2013-06-18 Thread Berend Hasselman
On 18-06-2013, at 13:23, peter dalgaard 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

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

2013-06-18 Thread Simone Giannerini
1. OpenSuse 12.1 R compiled against ACML 5.3.1 > sessionInfo() R version 3.0.1 RC (2013-05-08 r62723) Platform: x86_64-unknown-linux-gnu (64-bit) > A <- exp(-0.1*(row(jj)-col(jj))^2) > min(eigen(A,T,T)$values) [1] 2.521151e-10 > min(eigen(A+0i,T,T)$values) [1] 2.

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

2013-06-18 Thread Berend Hasselman
On 18-06-2013, at 09:57, peter dalgaard wrote: > > On Jun 18, 2013, at 03:30 , robin hankin wrote: > >> R-3.0.1 rev 62743, binary downloaded from CRAN just now; macosx 10.8.3 >> >> Hello, >> >> eigen(symmetric=TRUE) behaves strangely when given complex matrices. >> >> >> The following two

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

2013-06-18 Thread Davor Cubranic
On 13-06-18 12:57 AM, peter dalgaard wrote: Coercing A to a complex matrix should make no difference, but makes >eigen() return the wrong answer: > >>min(eigen(A+0i,T,T)$values) >[1] -0.359347 >> > >This is very, very wrong. Yep. I see this also on 10.6/7 (Snow Leopard, Lion) and 3.0.x, but

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

2013-06-18 Thread peter dalgaard
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. At smaller dimensions, a curious effect is visible: >

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

2013-06-18 Thread peter dalgaard
On Jun 18, 2013, at 03:30 , robin hankin wrote: > R-3.0.1 rev 62743, binary downloaded from CRAN just now; macosx 10.8.3 > > Hello, > > eigen(symmetric=TRUE) behaves strangely when given complex matrices. > > > The following two lines define 'A', a 100x100 (real) symmetric matrix > which theo

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

2013-06-17 Thread robin hankin
R-3.0.1 rev 62743, binary downloaded from CRAN just now; macosx 10.8.3 Hello, eigen(symmetric=TRUE) behaves strangely when given complex matrices. The following two lines define 'A', a 100x100 (real) symmetric matrix which theoretical considerations [Bochner's theorem] show to be positive defin