I just got a new twin processor Xeon machine and installed S-plus, only to find
that my calls to chol() don't work any more because t(x)%*%x is not a symmetric
matrix on the new machine.
Has anyone else seen this problem? Are there some hardware floating point
settings, or blas libraries that I need to tweak?
A similar thing happens with crossprod(), which could cause lm.chol() to fail.
Details follow, if anyone wants them.
Both machines running S-plus 6.1
# On old machine (Compaq Evo W4000 Pentium 4) (Windows 2000)
> getenv("PROCESSOR_IDENTIFIER")
PROCESSOR_IDENTIFIER
"x86 Family 15 Model 0 Stepping 10, GenuineIntel"
> getenv("PROCESSOR_LEVEL")
PROCESSOR_LEVEL
"15"
> getenv("PROCESSOR_REVISION")
PROCESSOR_REVISION
"000a"
> set.seed(1)
> x1 <- cbind(1,rnorm(100),rnorm(100))
> txx1 <- t(x1) %*% x1
> txx1 == t(txx1)
[,1] [,2] [,3]
[1,] T T T
[2,] T T T
[3,] T T T
> print(txx1, digits=17)
[,1] [,2] [,3]
[1,] 100.0000000000000000 -3.9481730102059434 -7.0954864056440785
[2,] -3.9481730102059434 92.2616446624295890 4.7666434140282590
[3,] -7.0954864056440785 4.7666434140282590 106.5686545795114500
>
# On new machine (HP xw8000 twin Xeon processor) (Windows 2000)
> getenv("PROCESSOR_IDENTIFIER")
PROCESSOR_IDENTIFIER
"x86 Family 15 Model 2 Stepping 7, GenuineIntel"
> getenv("PROCESSOR_LEVEL")
PROCESSOR_LEVEL
"15"
> getenv("PROCESSOR_REVISION")
PROCESSOR_REVISION
"0207"
> set.seed(1)
> x2 <- cbind(1,rnorm(100),rnorm(100))
> txx2 <- t(x2) %*% x2
> txx2 == t(txx2) ##### Should be symmetric, but is not !!!!!!!!
[,1] [,2] [,3]
[1,] T F F
[2,] F T T
[3,] F T T
> print(txx2, digits=20)
[,1] [,2] [,3]
[1,] 100.0000000000000000 -3.9481730102059425 -7.0954864056440794
[2,] -3.9481730102059434 92.2616446624295890 4.7666434140282590
[3,] -7.0954864056440785 4.7666434140282590 106.5686545795114500
> # Similar thing happens with crossprod()
> crossprod(x2, x2) == t(crossprod(x2,x2))
[,1] [,2] [,3]
[1,] T F T
[2,] F T T
[3,] T T T
>
# Do some comparisons of the objects
# The objects generated by rnorm() are identical on the two machines, but the
# results of matrix multiply are not
> all(x2==x1)
[1] T
> txx2==txx1
[,1] [,2] [,3]
[1,] T F F
[2,] T T T
[3,] T T T
> txx1 - txx2
[,1] [,2] [,3]
[1,] 0 -8.881784e-016 8.881784e-016
[2,] 0 0.000000e+000 0.000000e+000
[3,] 0 0.000000e+000 0.000000e+000
|