An embedded and charset-unspecified text was scrubbed... Name: not available URL: <https://stat.ethz.ch/pipermail/r-devel/attachments/20130618/af353d2a/attachment.pl>
eigen(symmetric=TRUE) for complex matrices
15 messages · robin hankin, Ravi Varadhan, Davor Cubranic +4 more
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 theoretical considerations [Bochner's theorem] show to be positive definite: jj <- matrix(0,100,100) A <- exp(-0.1*(row(jj)-col(jj))^2) A's being positive-definite is important to me:
min(eigen(A,T,T)$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,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 NOT with a MacPorts build of 2.15.3 that I had lying around. So this sits somewhere between Mac builds, R versions, and possibly LAPACK issues. Can anyone reproduce on non-Mac? It's not the first time complex matrices have acted up. We may need to put in a regression test to catch it earlier.
I would expect these two commands to return identical values, up to numerical precision. Compare svd():
dput(min(eigen(A,T,T)$values))
2.52115250343783e-10
dput(min(eigen(A+0i,T,T)$values))
-0.359346984206908
dput(min(svd(A)$d))
2.52115166468044e-10
dput(min(svd(A+0i)$d))
2.52115166468044e-10
So svd() doesn't care about the coercion to complex. The 'A' matrix isn't particularly badly conditioned:
eigen(A,T)$vectors -> e crossprod(e)[1:4,1:4]
also:
crossprod(A,solve(A))
[and the associated commands with A+0i in place of A], give errors of order 1e-14 or less. I think the eigenvectors are misbehaving too:
eigen(A,T)$vectors -> ev1 eigen(A+0i,T)$vectors -> ev2 range(Re((A %*% ev1[,100])/ev1[,100]))
[1] 2.497662e-10 2.566555e-10 # min=max mathematically; differences due to numerics
range(Re((A %*% ev2[,100])/ev2[,100]))
[1] -19.407290 4.412938 # off the scale errors [note the difference in sign]
FWIW, these problems do not appear to occur if symmetric=FALSE:
min(Re(eigen(A+0i,F,T)$values))
[1] 2.521153e-10
min(Re(eigen(A,F,T)$values))
[1] 2.521153e-10
and the eigenvectors appear to behave themselves too. Also, can I raise a doco? The documentation for eigen() is not entirely transparent with regard to the 'symmetric' argument. For complex matrices, 'symmetric' should read 'Hermitian':
B <- matrix(c(2,1i,-1i,2),2,2) # 'B' is Hermitian eigen(B,F,T)$values
[1] 3+0i 1+0i
eigen(B,T,T)$values # answers agree as expected if 'symmetric' means
'Hermitian' [1] 3 1
C <- matrix(c(2,1i,1i,2),2,2) # 'C' is symmetric eigen(C,F,T)$values
[1] 2-1i 2+1i
eigen(C,T,T)$values # answers disagree because 'C' is not Hermitian
[1] 3 1
This issue, however, I do not understand; ?eigen already has:
symmetric: if ?TRUE?, the matrix is assumed to be symmetric (or
Hermitian if complex) and only its lower triangle (diagonal
included) is used. If ?symmetric? is not specified, the
matrix is inspected for symmetry.
Which part of "Hermitian if complex" is unclear?
Peter Dalgaard, Professor, Center for Statistics, Copenhagen Business School Solbjerg Plads 3, 2000 Frederiksberg, Denmark Phone: (+45)38153501 Email: pd.mes at cbs.dk Priv: PDalgd at gmail.com
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:
jj <- matrix(0,33,33) A <- exp(-0.1*(row(jj)-col(jj))^2) eigen(A,,T)
$values [1] 5.497434e+00 5.187095e+00 4.708628e+00 4.112589e+00 3.456602e+00 [6] 2.796253e+00 2.177674e+00 1.633085e+00 1.179630e+00 8.209787e-01 [11] 5.506744e-01 3.560858e-01 2.220308e-01 1.335197e-01 7.744384e-02 [16] 4.332359e-02 2.337139e-02 1.215408e-02 6.089878e-03 2.937761e-03 [21] 1.363015e-03 6.073985e-04 2.595197e-04 1.060717e-04 4.134966e-05 [26] 1.531412e-05 5.360284e-06 1.760457e-06 5.369131e-07 1.496798e-07 [31] 3.715172e-08 7.807420e-09 1.238432e-09 $vectors NULL
eigen(A+0i,,T)
$values [1] 5.497434e+00 5.187095e+00 4.708628e+00 4.112589e+00 3.456602e+00 [6] 2.796253e+00 2.177674e+00 1.633085e+00 1.179630e+00 1.000000e+00 [11] 8.209787e-01 5.506744e-01 3.560858e-01 2.220308e-01 1.335197e-01 [16] 7.744384e-02 4.332359e-02 2.337139e-02 1.215408e-02 6.089878e-03 [21] 2.937761e-03 1.363015e-03 6.073985e-04 2.595197e-04 1.060717e-04 [26] 4.134966e-05 1.531412e-05 5.360284e-06 1.760457e-06 5.369131e-07 [31] 1.496793e-07 3.707403e-08 6.664029e-09 $vectors NULL Notice that a bogus eigenvalue of 1.000 has sneaked into the 10th position in the complex case. This seems to be happening with increasing frequency as the dimension is increased, e.g.
jj <- matrix(0,39,39) A <- exp(-0.1*(row(jj)-col(jj))^2) eigen(A+0i,,T)
$values [1] 5.525963e+00 5.295552e+00 4.932840e+00 4.466692e+00 4.313221e+00 [6] 3.931940e+00 3.365090e+00 2.800272e+00 2.266053e+00 1.992364e+00 [11] 1.783467e+00 1.365369e+00 1.016934e+00 7.369943e-01 5.730933e-01 [16] 5.197947e-01 3.568288e-01 2.384533e-01 1.551331e-01 1.070839e-01 [21] 9.826254e-02 6.059788e-02 3.638230e-02 2.126343e-02 1.318810e-02 [26] 1.209480e-02 6.693571e-03 3.602773e-03 1.884985e-03 1.011671e-03 [31] 9.580352e-04 4.725998e-04 2.260430e-04 1.046915e-04 4.687560e-05 [36] 3.794374e-05 2.020133e-05 7.751127e-06 1.472817e-06 $vectors NULL
eigen(A,,T)
$values [1] 5.525963e+00 5.295552e+00 4.932840e+00 4.466692e+00 3.931940e+00 [6] 3.365090e+00 2.800272e+00 2.266053e+00 1.783467e+00 1.365369e+00 [11] 1.016934e+00 7.369943e-01 5.197947e-01 3.568288e-01 2.384533e-01 [16] 1.551331e-01 9.826254e-02 6.059788e-02 3.638230e-02 2.126343e-02 [21] 1.209480e-02 6.693571e-03 3.602773e-03 1.884985e-03 9.580352e-04 [26] 4.725998e-04 2.260430e-04 1.046915e-04 4.687623e-05 2.025048e-05 [31] 8.418654e-06 3.356873e-06 1.278239e-06 4.620504e-07 1.572170e-07 [36] 4.971980e-08 1.431462e-08 3.613385e-09 7.510902e-10 $vectors NULL in which most of the rightmost column of the complex case appear to be insertions. -pd
On Jun 18, 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 lines define 'A', a 100x100 (real) symmetric matrix which theoretical considerations [Bochner's theorem] show to be positive definite: jj <- matrix(0,100,100) A <- exp(-0.1*(row(jj)-col(jj))^2) A's being positive-definite is important to me:
min(eigen(A,T,T)$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,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 NOT with a MacPorts build of 2.15.3 that I had lying around. So this sits somewhere between Mac builds, R versions, and possibly LAPACK issues. Can anyone reproduce on non-Mac? It's not the first time complex matrices have acted up. We may need to put in a regression test to catch it earlier.
I would expect these two commands to return identical values, up to numerical precision. Compare svd():
dput(min(eigen(A,T,T)$values))
2.52115250343783e-10
dput(min(eigen(A+0i,T,T)$values))
-0.359346984206908
dput(min(svd(A)$d))
2.52115166468044e-10
dput(min(svd(A+0i)$d))
2.52115166468044e-10
So svd() doesn't care about the coercion to complex. The 'A' matrix isn't particularly badly conditioned:
eigen(A,T)$vectors -> e crossprod(e)[1:4,1:4]
also:
crossprod(A,solve(A))
[and the associated commands with A+0i in place of A], give errors of order 1e-14 or less. I think the eigenvectors are misbehaving too:
eigen(A,T)$vectors -> ev1 eigen(A+0i,T)$vectors -> ev2 range(Re((A %*% ev1[,100])/ev1[,100]))
[1] 2.497662e-10 2.566555e-10 # min=max mathematically; differences due to numerics
range(Re((A %*% ev2[,100])/ev2[,100]))
[1] -19.407290 4.412938 # off the scale errors [note the difference in sign]
FWIW, these problems do not appear to occur if symmetric=FALSE:
min(Re(eigen(A+0i,F,T)$values))
[1] 2.521153e-10
min(Re(eigen(A,F,T)$values))
[1] 2.521153e-10
and the eigenvectors appear to behave themselves too. Also, can I raise a doco? The documentation for eigen() is not entirely transparent with regard to the 'symmetric' argument. For complex matrices, 'symmetric' should read 'Hermitian':
B <- matrix(c(2,1i,-1i,2),2,2) # 'B' is Hermitian eigen(B,F,T)$values
[1] 3+0i 1+0i
eigen(B,T,T)$values # answers agree as expected if 'symmetric' means
'Hermitian' [1] 3 1
C <- matrix(c(2,1i,1i,2),2,2) # 'C' is symmetric eigen(C,F,T)$values
[1] 2-1i 2+1i
eigen(C,T,T)$values # answers disagree because 'C' is not Hermitian
[1] 3 1
This issue, however, I do not understand; ?eigen already has:
symmetric: if ?TRUE?, the matrix is assumed to be symmetric (or
Hermitian if complex) and only its lower triangle (diagonal
included) is used. If ?symmetric? is not specified, the
matrix is inspected for symmetry.
Which part of "Hermitian if complex" is unclear?
--
Peter Dalgaard, Professor,
Center for Statistics, Copenhagen Business School
Solbjerg Plads 3, 2000 Frederiksberg, Denmark
Phone: (+45)38153501
Email: pd.mes at cbs.dk Priv: PDalgd at gmail.com
Peter Dalgaard, Professor Center for Statistics, Copenhagen Business School Solbjerg Plads 3, 2000 Frederiksberg, Denmark Phone: (+45)38153501 Email: pd.mes at cbs.dk Priv: PDalgd at gmail.com
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
version
_ platform i386-w64-mingw32 arch i386 os mingw32 system i386, mingw32 status major 3 minor 0.1 year 2013 month 05 day 16 svn rev 62743 language R version.string R version 3.0.1 (2013-05-16) nickname Good Sport
-----Original Message----- From: r-devel-bounces at r-project.org [mailto:r-devel-bounces at r-project.org] On Behalf Of peter dalgaard Sent: Tuesday, June 18, 2013 7:23 AM To: robin hankin Cc: r-devel at r-project.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 33x33 and above. At smaller dimensions, a curious effect is visible:
jj <- matrix(0,33,33) A <- exp(-0.1*(row(jj)-col(jj))^2) eigen(A,,T)
$values [1] 5.497434e+00 5.187095e+00 4.708628e+00 4.112589e+00 3.456602e+00 [6] 2.796253e+00 2.177674e+00 1.633085e+00 1.179630e+00 8.209787e-01 [11] 5.506744e-01 3.560858e-01 2.220308e-01 1.335197e-01 7.744384e-02 [16] 4.332359e-02 2.337139e-02 1.215408e-02 6.089878e-03 2.937761e-03 [21] 1.363015e-03 6.073985e-04 2.595197e-04 1.060717e-04 4.134966e-05 [26] 1.531412e-05 5.360284e-06 1.760457e-06 5.369131e-07 1.496798e-07 [31] 3.715172e-08 7.807420e-09 1.238432e-09 $vectors NULL
eigen(A+0i,,T)
$values [1] 5.497434e+00 5.187095e+00 4.708628e+00 4.112589e+00 3.456602e+00 [6] 2.796253e+00 2.177674e+00 1.633085e+00 1.179630e+00 1.000000e+00 [11] 8.209787e-01 5.506744e-01 3.560858e-01 2.220308e-01 1.335197e-01 [16] 7.744384e-02 4.332359e-02 2.337139e-02 1.215408e-02 6.089878e-03 [21] 2.937761e-03 1.363015e-03 6.073985e-04 2.595197e-04 1.060717e-04 [26] 4.134966e-05 1.531412e-05 5.360284e-06 1.760457e-06 5.369131e-07 [31] 1.496793e-07 3.707403e-08 6.664029e-09 $vectors NULL Notice that a bogus eigenvalue of 1.000 has sneaked into the 10th position in the complex case. This seems to be happening with increasing frequency as the dimension is increased, e.g.
jj <- matrix(0,39,39) A <- exp(-0.1*(row(jj)-col(jj))^2) eigen(A+0i,,T)
$values [1] 5.525963e+00 5.295552e+00 4.932840e+00 4.466692e+00 4.313221e+00 [6] 3.931940e+00 3.365090e+00 2.800272e+00 2.266053e+00 1.992364e+00 [11] 1.783467e+00 1.365369e+00 1.016934e+00 7.369943e-01 5.730933e-01 [16] 5.197947e-01 3.568288e-01 2.384533e-01 1.551331e-01 1.070839e-01 [21] 9.826254e-02 6.059788e-02 3.638230e-02 2.126343e-02 1.318810e-02 [26] 1.209480e-02 6.693571e-03 3.602773e-03 1.884985e-03 1.011671e-03 [31] 9.580352e-04 4.725998e-04 2.260430e-04 1.046915e-04 4.687560e-05 [36] 3.794374e-05 2.020133e-05 7.751127e-06 1.472817e-06 $vectors NULL
eigen(A,,T)
$values [1] 5.525963e+00 5.295552e+00 4.932840e+00 4.466692e+00 3.931940e+00 [6] 3.365090e+00 2.800272e+00 2.266053e+00 1.783467e+00 1.365369e+00 [11] 1.016934e+00 7.369943e-01 5.197947e-01 3.568288e-01 2.384533e-01 [16] 1.551331e-01 9.826254e-02 6.059788e-02 3.638230e-02 2.126343e-02 [21] 1.209480e-02 6.693571e-03 3.602773e-03 1.884985e-03 9.580352e-04 [26] 4.725998e-04 2.260430e-04 1.046915e-04 4.687623e-05 2.025048e-05 [31] 8.418654e-06 3.356873e-06 1.278239e-06 4.620504e-07 1.572170e-07 [36] 4.971980e-08 1.431462e-08 3.613385e-09 7.510902e-10 $vectors NULL in which most of the rightmost column of the complex case appear to be insertions. -pd
On Jun 18, 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 lines define 'A', a 100x100 (real) symmetric matrix which theoretical considerations [Bochner's theorem] show to be positive definite: jj <- matrix(0,100,100) A <- exp(-0.1*(row(jj)-col(jj))^2) A's being positive-definite is important to me:
min(eigen(A,T,T)$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,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 NOT with a MacPorts build of 2.15.3 that I had lying around. So this sits somewhere between Mac builds, R versions, and possibly LAPACK issues. Can anyone reproduce on non-Mac? It's not the first time complex matrices have acted up. We may need to put in a regression test to catch it earlier.
I would expect these two commands to return identical values, up to numerical precision. Compare svd():
dput(min(eigen(A,T,T)$values))
2.52115250343783e-10
dput(min(eigen(A+0i,T,T)$values))
-0.359346984206908
dput(min(svd(A)$d))
2.52115166468044e-10
dput(min(svd(A+0i)$d))
2.52115166468044e-10
So svd() doesn't care about the coercion to complex. The 'A' matrix isn't particularly badly conditioned:
eigen(A,T)$vectors -> e crossprod(e)[1:4,1:4]
also:
crossprod(A,solve(A))
[and the associated commands with A+0i in place of A], give errors of order 1e-14 or less. I think the eigenvectors are misbehaving too:
eigen(A,T)$vectors -> ev1 eigen(A+0i,T)$vectors -> ev2 range(Re((A %*% ev1[,100])/ev1[,100]))
[1] 2.497662e-10 2.566555e-10 # min=max mathematically; differences due to numerics
range(Re((A %*% ev2[,100])/ev2[,100]))
[1] -19.407290 4.412938 # off the scale errors [note the difference in sign]
FWIW, these problems do not appear to occur if symmetric=FALSE:
min(Re(eigen(A+0i,F,T)$values))
[1] 2.521153e-10
min(Re(eigen(A,F,T)$values))
[1] 2.521153e-10
and the eigenvectors appear to behave themselves too. Also, can I raise a doco? The documentation for eigen() is not entirely transparent with regard to the 'symmetric' argument. For complex matrices, 'symmetric' should read 'Hermitian':
B <- matrix(c(2,1i,-1i,2),2,2) # 'B' is Hermitian eigen(B,F,T)$values
[1] 3+0i 1+0i
eigen(B,T,T)$values # answers agree as expected if 'symmetric' means
'Hermitian' [1] 3 1
C <- matrix(c(2,1i,1i,2),2,2) # 'C' is symmetric eigen(C,F,T)$values
[1] 2-1i 2+1i
eigen(C,T,T)$values # answers disagree because 'C' is not Hermitian
[1] 3 1
This issue, however, I do not understand; ?eigen already has:
symmetric: if 'TRUE', the matrix is assumed to be symmetric (or
Hermitian if complex) and only its lower triangle (diagonal
included) is used. If 'symmetric' is not specified, the
matrix is inspected for symmetry.
Which part of "Hermitian if complex" is unclear?
--
Peter Dalgaard, Professor,
Center for Statistics, Copenhagen Business School Solbjerg Plads 3,
2000 Frederiksberg, Denmark
Phone: (+45)38153501
Email: pd.mes at cbs.dk Priv: PDalgd at gmail.com
-- Peter Dalgaard, Professor Center for Statistics, Copenhagen Business School Solbjerg Plads 3, 2000 Frederiksberg, Denmark Phone: (+45)38153501 Email: pd.mes at cbs.dk Priv: PDalgd at gmail.com ______________________________________________ R-devel at r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-devel
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 NOT with a MacPorts build of 2.15.3 that I had lying around. So this sits somewhere between Mac builds, R versions, and possibly LAPACK issues. Can anyone reproduce on non-Mac?
It doesn't happen with the CRAN binary of 2.15.3 built for Ubuntu Precise. Davor
On 18-06-2013, at 09:57, peter dalgaard <pdalgd at gmail.com> 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 lines define 'A', a 100x100 (real) symmetric matrix which theoretical considerations [Bochner's theorem] show to be positive definite: jj <- matrix(0,100,100) A <- exp(-0.1*(row(jj)-col(jj))^2) A's being positive-definite is important to me:
min(eigen(A,T,T)$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,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 NOT with a MacPorts build of 2.15.3 that I had lying around. So this sits somewhere between Mac builds, R versions, and possibly LAPACK issues. Can anyone reproduce on non-Mac?
The problem does not occur with the Cran binary of R-3.0.1 Kubuntu 12.04 64-bit. That R uses the system provided Blas (libblas 1.2.30110419) and Lapack 3.3.1; I don't know if these have been patched. I have been able to reproduce the problem on a self compiled version of R-3.0.1 using Rlapack and Rblas on Ubuntu 10.04 64-bit. Berend
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.521154e-10
2. OpenSuse 12.3 R compiled against Intel MKL 10.2
R version 3.0.0 Patched (2013-04-23 r62650)
Platform: x86_64-unknown-linux-gnu (64-bit)
> jj <- matrix(0,100,100)
> A <- exp(-0.1*(row(jj)-col(jj))^2)
> min(eigen(A,T,T)$values)
[1] 2.521153e-10
> min(eigen(A+0i,T,T)$values)
[1] 2.521154e-10
3. Windows 7 64, R 3.0.0p (binary, built-in libraries)
R version 3.0.0 Patched (2013-04-03 r62488)
Platform: x86_64-w64-mingw32/x64 (64-bit)
> jj <- matrix(0,100,100)
> A <- exp(-0.1*(row(jj)-col(jj))^2)
> min(eigen(A,T,T)$values)
[1] 2.521153e-10
> min(eigen(A+0i,T,T)$values)
[1] -0.359347
same behaviour with the 32 bit binaries.
On Tue, Jun 18, 2013 at 9:57 AM, peter dalgaard <pdalgd at gmail.com> 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 lines define 'A', a 100x100 (real) symmetric matrix which theoretical considerations [Bochner's theorem] show to be positive definite: jj <- matrix(0,100,100) A <- exp(-0.1*(row(jj)-col(jj))^2) A's being positive-definite is important to me:
min(eigen(A,T,T)$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,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 NOT with a MacPorts build of 2.15.3 that I had lying around. So this sits somewhere between Mac builds, R versions, and possibly LAPACK issues. Can anyone reproduce on non-Mac? It's not the first time complex matrices have acted up. We may need to put in a regression test to catch it earlier.
I would expect these two commands to return identical values, up to numerical precision. Compare svd():
dput(min(eigen(A,T,T)$values))
2.52115250343783e-10
dput(min(eigen(A+0i,T,T)$values))
-0.359346984206908
dput(min(svd(A)$d))
2.52115166468044e-10
dput(min(svd(A+0i)$d))
2.52115166468044e-10
So svd() doesn't care about the coercion to complex. The 'A' matrix isn't particularly badly conditioned:
eigen(A,T)$vectors -> e crossprod(e)[1:4,1:4]
also:
crossprod(A,solve(A))
[and the associated commands with A+0i in place of A], give errors of order 1e-14 or less. I think the eigenvectors are misbehaving too:
eigen(A,T)$vectors -> ev1 eigen(A+0i,T)$vectors -> ev2 range(Re((A %*% ev1[,100])/ev1[,100]))
[1] 2.497662e-10 2.566555e-10 # min=max mathematically; differences due to numerics
range(Re((A %*% ev2[,100])/ev2[,100]))
[1] -19.407290 4.412938 # off the scale errors [note the difference in sign]
FWIW, these problems do not appear to occur if symmetric=FALSE:
min(Re(eigen(A+0i,F,T)$values))
[1] 2.521153e-10
min(Re(eigen(A,F,T)$values))
[1] 2.521153e-10
and the eigenvectors appear to behave themselves too. Also, can I raise a doco? The documentation for eigen() is not entirely transparent with regard to the 'symmetric' argument. For complex matrices, 'symmetric' should read 'Hermitian':
B <- matrix(c(2,1i,-1i,2),2,2) # 'B' is Hermitian eigen(B,F,T)$values
[1] 3+0i 1+0i
eigen(B,T,T)$values # answers agree as expected if 'symmetric' means
'Hermitian' [1] 3 1
C <- matrix(c(2,1i,1i,2),2,2) # 'C' is symmetric eigen(C,F,T)$values
[1] 2-1i 2+1i
eigen(C,T,T)$values # answers disagree because 'C' is not Hermitian
[1] 3 1
This issue, however, I do not understand; ?eigen already has:
symmetric: if ?TRUE?, the matrix is assumed to be symmetric (or
Hermitian if complex) and only its lower triangle (diagonal
included) is used. If ?symmetric? is not specified, the
matrix is inspected for symmetry.
Which part of "Hermitian if complex" is unclear?
--
Peter Dalgaard, Professor,
Center for Statistics, Copenhagen Business School
Solbjerg Plads 3, 2000 Frederiksberg, Denmark
Phone: (+45)38153501
Email: pd.mes at cbs.dk Priv: PDalgd at gmail.com
______________________________________________ R-devel at r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-devel
-- ___________________________________________________ Simone Giannerini Dipartimento di Scienze Statistiche "Paolo Fortunati" Universita' di Bologna Via delle belle arti 41 - 40126 Bologna, ITALY Tel: +39 051 2098262 Fax: +39 051 232153 http://www2.stat.unibo.it/giannerini/
On 18-06-2013, at 13:23, peter dalgaard <pdalgd at 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
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 apparent pattern is that things work when R is linked against external libraries and not when the internal ones are used. So it could be time to start looking for differences between R's lapack module and the original LAPACK code. -pd
Peter Dalgaard, Professor, Center for Statistics, Copenhagen Business School Solbjerg Plads 3, 2000 Frederiksberg, Denmark Phone: (+45)38153501 Email: pd.mes at cbs.dk Priv: PDalgd at gmail.com
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 finger seems to be pointing to the internal libRlapack/libRblas. The apparent pattern is that things work when R is linked against external libraries and not when the internal ones are used. So it could be time to start looking for differences between R's lapack module and the original LAPACK code.
Now done and no (relevant) differences found. Hmmm.... 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. One option could be to bypass R by coding Robin's example in pure Fortran and see whether issues persist when linked against libRlapack, vecLib, and the relevant subset of current LAPACK sources. (The bogus 1.0000 eigenvalue in the 33x33 variant of Robin's example should make it rather easy to spot whether things work or not.)
Peter Dalgaard, Professor, Center for Statistics, Copenhagen Business School Solbjerg Plads 3, 2000 Frederiksberg, Denmark Phone: (+45)38153501 Email: pd.mes at cbs.dk Priv: PDalgd at gmail.com
On 19-06-2013, at 10:24, peter dalgaard <pdalgd at gmail.com> 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 algorithm gives the correct result.
Thanks Berend, The finger seems to be pointing to the internal libRlapack/libRblas. The apparent pattern is that things work when R is linked against external libraries and not when the internal ones are used. So it could be time to start looking for differences between R's lapack module and the original LAPACK code.
Now done and no (relevant) differences found. Hmmm.... 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. One option could be to bypass R by coding Robin's example in pure Fortran and see whether issues persist when linked against libRlapack, vecLib, and the relevant subset of current LAPACK sources. (The bogus 1.0000 eigenvalue in the 33x33 variant of Robin's example should make it rather easy to spot whether things work or not.)
I have made that pure Fortran program.
I have run it on OS X 10.8.4 with
- libRblas libRlapack
- my private refblas and Lapack
- framework Accelerate
- downloaded zheev.f + dependencies + libRblas
- downloaded zheev.f + dependencies + my refblas
The Fortran program and the shell script to do the compiling and running and the output file follow at the end of this mail.
The shell scripts runs otool -L on each executable to make sure it is linked to the intended libraries.
The fortran program runs zheev with the minimal lwork and the optimal lwork.
Summary of the output:
- libRblas libRlapack ==> Bad
- my private refblas and Lapack ==> OK
- framework Accelerate ==> OK
- downloaded zheev.f + dependencies + libRblas ==> Bad
- downloaded zheev.f + dependencies + my refblas ==> OK
It looks like libRblas is the cause of the problem.
I haven't done any further investigation of the matter.
Berend
Output is:
<output>
Running hb with Rlapack and Rblas
./hbRlapack:
/Library/Frameworks/R.framework/Versions/3.0/Resources/lib/libRblas.dylib (compatibility version 0.0.0, current version 0.0.0)
/Library/Frameworks/R.framework/Versions/3.0/Resources/lib/libRlapack.dylib (compatibility version 3.0.0, current version 3.0.1)
/usr/lib/libSystem.B.dylib (compatibility version 1.0.0, current version 169.3.0)
lwork= 199
info= 0
min eig value= -4.433595556845234E-008
lwork= 3300
info= 0
min eig value= -0.359347003948635
Running hb with Haslapack and Hasblas (Lapack 3.4.2)
./hbHaslapack:
/Users/berendhasselman/Library/lib/x86_64/librefblas.dylib (compatibility version 0.0.0, current version 0.0.0)
/Users/berendhasselman/Library/lib/x86_64/liblapack.dylib (compatibility version 0.0.0, current version 0.0.0)
/usr/lib/libSystem.B.dylib (compatibility version 1.0.0, current version 169.3.0)
lwork= 199
info= 0
min eig value= -4.433595556845234E-008
lwork= 3300
info= 0
min eig value= -4.433595559424842E-008
Running hb with -framework Accelerate
./hbveclib:
/System/Library/Frameworks/Accelerate.framework/Versions/A/Accelerate (compatibility version 1.0.0, current version 4.0.0)
/usr/lib/libSystem.B.dylib (compatibility version 1.0.0, current version 169.3.0)
lwork= 199
info= 0
min eig value= -4.433595568797253E-008
lwork= 3300
info= 0
min eig value= -4.433595550683581E-008
Compile downloaded zheev + dependencies
/Users/berendhasselman/Documents/Programming/R/problems/bug-error/eigbug
Running hb with downloaded zheev and Rblas
./hbzheev:
/Library/Frameworks/R.framework/Versions/3.0/Resources/lib/libRblas.dylib (compatibility version 0.0.0, current version 0.0.0)
/usr/lib/libSystem.B.dylib (compatibility version 1.0.0, current version 169.3.0)
lwork= 199
info= 0
min eig value= -4.433595556845234E-008
lwork= 3300
info= 0
min eig value= -0.359347003948635
Running hb with downloaded zheev and Hasrefblas
./hbzheev:
/Users/berendhasselman/Library/lib/x86_64/librefblas.dylib (compatibility version 0.0.0, current version 0.0.0)
/usr/lib/libSystem.B.dylib (compatibility version 1.0.0, current version 169.3.0)
lwork= 199
info= 0
min eig value= -4.433595556845234E-008
lwork= 3300
info= 0
min eig value= -4.433595559424842E-008
</output>
Fortran program:
<fortran>
program test
integer n
parameter(n=100)
complex*16 A(n,n)
double precision w(n), rwork(3*n-2)
integer lwork, info
complex*16 work, wtmp(1)
allocatable work(:)
call genmat(A,n)
lwork = 2*n-1
allocate(work(lwork))
print *, "lwork=",lwork
call zheev("N","L",n,A,n,w,work,lwork,rwork,info)
print *, "info=",info
print *, "min eig value=",minval(w)
deallocate(work)
call genmat(A,n)
lwork = -1
call zheev("N","L",n,A,n,w,wtmp,lwork,rwork,info)
lwork = real(wtmp(1))
allocate(work(lwork))
print *, "lwork=",lwork
call zheev("N","L",n,A,n,w,work,lwork,rwork,info)
print *, "info=",info
print *, "min eig value=",minval(w)
stop
end
subroutine genmat(A,n)
integer n
complex*16 A(n,*)
integer i, j
double precision tmp
do i=1,n
do j=1,n
tmp = exp(-0.1D0 * (i-j)**2)
A(i,j) = cmplx(tmp,0.0D0)
enddo
enddo
return
end
</fortran>
Shell script to do the compiling and running:
<script>
gfortran -c hankinbug.f -o hankinbug.o
gfortran hankinbug.o -L/Library/Frameworks/R.framework/Libraries -lRblas -lRlapack -o hbRlapack
echo "Running hb with Rlapack and Rblas"
otool -L ./hbRlapack
./hbRlapack
echo ""
gfortran hankinbug.o -L${HOME}/Library/lib/x86_64 -lrefblas -llapack -o hbHaslapack
echo "Running hb with Haslapack and Hasblas (Lapack 3.4.2)"
otool -L ./hbHaslapack
./hbHaslapack
echo ""
gfortran hankinbug.o -framework Accelerate -o hbveclib
echo "Running hb with -framework Accelerate"
otool -L ./hbveclib
./hbveclib
echo ""
echo "Compile downloaded zheev + dependencies"
cd zheev
gfortran -c *.f
cd -
gfortran hankinbug.o zheev/*.o -L/Library/Frameworks/R.framework/Libraries -lRblas -o hbzheev
echo "Running hb with downloaded zheev and Rblas"
otool -L ./hbzheev
./hbzheev
echo ""
gfortran hankinbug.o zheev/*.o -L${HOME}/Library/lib/x86_64 -lrefblas -o hbzheev
echo "Running hb with downloaded zheev and Hasrefblas"
otool -L ./hbzheev
./hbzheev
echo ""
</script>
On Jun 19, 2013, at 12:43 , Berend Hasselman wrote:
On 19-06-2013, at 10:24, peter dalgaard <pdalgd at gmail.com> 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 algorithm gives the correct result.
Thanks Berend, The finger seems to be pointing to the internal libRlapack/libRblas. The apparent pattern is that things work when R is linked against external libraries and not when the internal ones are used. So it could be time to start looking for differences between R's lapack module and the original LAPACK code.
Now done and no (relevant) differences found. Hmmm.... 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. One option could be to bypass R by coding Robin's example in pure Fortran and see whether issues persist when linked against libRlapack, vecLib, and the relevant subset of current LAPACK sources. (The bogus 1.0000 eigenvalue in the 33x33 variant of Robin's example should make it rather easy to spot whether things work or not.)
I have made that pure Fortran program.
I have run it on OS X 10.8.4 with
- libRblas libRlapack
- my private refblas and Lapack
- framework Accelerate
- downloaded zheev.f + dependencies + libRblas
- downloaded zheev.f + dependencies + my refblas
The Fortran program and the shell script to do the compiling and running and the output file follow at the end of this mail.
The shell scripts runs otool -L on each executable to make sure it is linked to the intended libraries.
The fortran program runs zheev with the minimal lwork and the optimal lwork.
Summary of the output:
- libRblas libRlapack ==> Bad
- my private refblas and Lapack ==> OK
- framework Accelerate ==> OK
- downloaded zheev.f + dependencies + libRblas ==> Bad
- downloaded zheev.f + dependencies + my refblas ==> OK
It looks like libRblas is the cause of the problem.
I haven't done any further investigation of the matter.
Berend
Output is:
<output>
Running hb with Rlapack and Rblas
./hbRlapack:
/Library/Frameworks/R.framework/Versions/3.0/Resources/lib/libRblas.dylib (compatibility version 0.0.0, current version 0.0.0)
/Library/Frameworks/R.framework/Versions/3.0/Resources/lib/libRlapack.dylib (compatibility version 3.0.0, current version 3.0.1)
/usr/lib/libSystem.B.dylib (compatibility version 1.0.0, current version 169.3.0)
lwork= 199
info= 0
min eig value= -4.433595556845234E-008
lwork= 3300
info= 0
min eig value= -0.359347003948635
Running hb with Haslapack and Hasblas (Lapack 3.4.2)
./hbHaslapack:
/Users/berendhasselman/Library/lib/x86_64/librefblas.dylib (compatibility version 0.0.0, current version 0.0.0)
/Users/berendhasselman/Library/lib/x86_64/liblapack.dylib (compatibility version 0.0.0, current version 0.0.0)
/usr/lib/libSystem.B.dylib (compatibility version 1.0.0, current version 169.3.0)
lwork= 199
info= 0
min eig value= -4.433595556845234E-008
lwork= 3300
info= 0
min eig value= -4.433595559424842E-008
Running hb with -framework Accelerate
./hbveclib:
/System/Library/Frameworks/Accelerate.framework/Versions/A/Accelerate (compatibility version 1.0.0, current version 4.0.0)
/usr/lib/libSystem.B.dylib (compatibility version 1.0.0, current version 169.3.0)
lwork= 199
info= 0
min eig value= -4.433595568797253E-008
lwork= 3300
info= 0
min eig value= -4.433595550683581E-008
Compile downloaded zheev + dependencies
/Users/berendhasselman/Documents/Programming/R/problems/bug-error/eigbug
Running hb with downloaded zheev and Rblas
./hbzheev:
/Library/Frameworks/R.framework/Versions/3.0/Resources/lib/libRblas.dylib (compatibility version 0.0.0, current version 0.0.0)
/usr/lib/libSystem.B.dylib (compatibility version 1.0.0, current version 169.3.0)
lwork= 199
info= 0
min eig value= -4.433595556845234E-008
lwork= 3300
info= 0
min eig value= -0.359347003948635
Running hb with downloaded zheev and Hasrefblas
./hbzheev:
/Users/berendhasselman/Library/lib/x86_64/librefblas.dylib (compatibility version 0.0.0, current version 0.0.0)
/usr/lib/libSystem.B.dylib (compatibility version 1.0.0, current version 169.3.0)
lwork= 199
info= 0
min eig value= -4.433595556845234E-008
lwork= 3300
info= 0
min eig value= -4.433595559424842E-008
</output>
Fortran program:
<fortran>
program test
integer n
parameter(n=100)
complex*16 A(n,n)
double precision w(n), rwork(3*n-2)
integer lwork, info
complex*16 work, wtmp(1)
allocatable work(:)
call genmat(A,n)
lwork = 2*n-1
allocate(work(lwork))
print *, "lwork=",lwork
call zheev("N","L",n,A,n,w,work,lwork,rwork,info)
print *, "info=",info
print *, "min eig value=",minval(w)
deallocate(work)
call genmat(A,n)
lwork = -1
call zheev("N","L",n,A,n,w,wtmp,lwork,rwork,info)
lwork = real(wtmp(1))
allocate(work(lwork))
print *, "lwork=",lwork
call zheev("N","L",n,A,n,w,work,lwork,rwork,info)
print *, "info=",info
print *, "min eig value=",minval(w)
stop
end
subroutine genmat(A,n)
integer n
complex*16 A(n,*)
integer i, j
double precision tmp
do i=1,n
do j=1,n
tmp = exp(-0.1D0 * (i-j)**2)
A(i,j) = cmplx(tmp,0.0D0)
enddo
enddo
return
end
</fortran>
Shell script to do the compiling and running:
<script>
gfortran -c hankinbug.f -o hankinbug.o
gfortran hankinbug.o -L/Library/Frameworks/R.framework/Libraries -lRblas -lRlapack -o hbRlapack
echo "Running hb with Rlapack and Rblas"
otool -L ./hbRlapack
./hbRlapack
echo ""
gfortran hankinbug.o -L${HOME}/Library/lib/x86_64 -lrefblas -llapack -o hbHaslapack
echo "Running hb with Haslapack and Hasblas (Lapack 3.4.2)"
otool -L ./hbHaslapack
./hbHaslapack
echo ""
gfortran hankinbug.o -framework Accelerate -o hbveclib
echo "Running hb with -framework Accelerate"
otool -L ./hbveclib
./hbveclib
echo ""
echo "Compile downloaded zheev + dependencies"
cd zheev
gfortran -c *.f
cd -
gfortran hankinbug.o zheev/*.o -L/Library/Frameworks/R.framework/Libraries -lRblas -o hbzheev
echo "Running hb with downloaded zheev and Rblas"
otool -L ./hbzheev
./hbzheev
echo ""
gfortran hankinbug.o zheev/*.o -L${HOME}/Library/lib/x86_64 -lrefblas -o hbzheev
echo "Running hb with downloaded zheev and Hasrefblas"
otool -L ./hbzheev
./hbzheev
echo ""
</script>
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 of the changes looked like this:
@@ -1561,7 +1561,7 @@
C( J, J ) = DBLE( C( J, J ) )
END IF
DO 170 L = 1, K
- IF( ( A( J, L ).NE.ZERO ) .OR. ( B( J, L ).NE.ZERO ) )
+c IF( ( A( J, L ).NE.ZERO ) .OR. ( B( J, L ).NE.ZERO ) )
$ THEN
TEMP1 = ALPHA*DCONJG( B( J, L ) )
TEMP2 = DCONJG( ALPHA*A( J, L ) )
Notice that the continuation line was NOT commented out. So FORTRAN, being what it is, continues the line before the comment and parses it as
DO 170 L = 1, KTHEN
with KTHEN uninitialized! and things go downhill from there.
(The uninitialized variable was actually hinted at in PR14964 and the fact that I could get one of my builds to segfault also helped.)
Peter Dalgaard, Professor Center for Statistics, Copenhagen Business School Solbjerg Plads 3, 2000 Frederiksberg, Denmark Phone: (+45)38153501 Email: pd.mes at cbs.dk Priv: PDalgd at gmail.com
peter dalgaard <pdalgd <at> 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,only.values=TRUE,symmetric=TRUE)$values) min(eigen(A+0i,only.values=TRUE,symmetric=TRUE)$values) with R's BLAS in 3.0.1 on Fedora 18 with gcc 4.7.2-8 gives the same error. Both 3.0.1 and R-devel (R's BLAS) on RHEL 5 gcc 4.1.2 give a segfault at zher2k(), in the 3.0.1 case at line 1566 in cmplxblas.f, at line 1570 for R-devel. So yes, there is a compiler issue too. The backtrace is La_rs_cmplx at Lapack.c:829 -> zheev at cmplx.f:10356 -> zhetrd at cmplx.f:11080 -> zher2k at cmplxblas.f:1570 for R-devel, and La_rs_cmplx at Lapack.c:829 -> zheev at cmplx.f:10356 -> zhetrd at cmplx.f:11080 -> zher2k at cmplxblas.f:1566 for 3.0.1. Roger
On 19-06-2013, at 14:17, peter dalgaard <pdalgd at gmail.com> 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 of the changes looked like this:
@@ -1561,7 +1561,7 @@
C( J, J ) = DBLE( C( J, J ) )
END IF
DO 170 L = 1, K
- IF( ( A( J, L ).NE.ZERO ) .OR. ( B( J, L ).NE.ZERO ) )
+c IF( ( A( J, L ).NE.ZERO ) .OR. ( B( J, L ).NE.ZERO ) )
$ THEN
TEMP1 = ALPHA*DCONJG( B( J, L ) )
TEMP2 = DCONJG( ALPHA*A( J, L ) )
Notice that the continuation line was NOT commented out. So FORTRAN, being what it is, continues the line before the comment and parses it as
DO 170 L = 1, KTHEN
with KTHEN uninitialized! and things go downhill from there.
(The uninitialized variable was actually hinted at in PR14964 and the fact that I could get one of my builds to segfault also helped.)
And it seems that line 1536 and 1537 ( R-3.0.1 source) have a similar change, so that will also have be put right. Berend
On Jun 19, 2013, at 15:58 , Berend Hasselman wrote:
On 19-06-2013, at 14:17, peter dalgaard <pdalgd at gmail.com> 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 of the changes looked like this:
@@ -1561,7 +1561,7 @@
C( J, J ) = DBLE( C( J, J ) )
END IF
DO 170 L = 1, K
- IF( ( A( J, L ).NE.ZERO ) .OR. ( B( J, L ).NE.ZERO ) )
+c IF( ( A( J, L ).NE.ZERO ) .OR. ( B( J, L ).NE.ZERO ) )
$ THEN
TEMP1 = ALPHA*DCONJG( B( J, L ) )
TEMP2 = DCONJG( ALPHA*A( J, L ) )
Notice that the continuation line was NOT commented out. So FORTRAN, being what it is, continues the line before the comment and parses it as
DO 170 L = 1, KTHEN
with KTHEN uninitialized! and things go downhill from there.
(The uninitialized variable was actually hinted at in PR14964 and the fact that I could get one of my builds to segfault also helped.)
And it seems that line 1536 and 1537 ( R-3.0.1 source) have a similar change, so that will also have be put right. Berend
Yes, that's the other one of the two changes. -pd
Peter Dalgaard, Professor Center for Statistics, Copenhagen Business School Solbjerg Plads 3, 2000 Frederiksberg, Denmark Phone: (+45)38153501 Email: pd.mes at cbs.dk Priv: PDalgd at gmail.com