Full_Name: Hao Yu
Version: 1.4.1
OS: Windws and Linux
Submission from: (NULL) (129.100.45.161)
Hi,
Recently we did some benchmarks regarding matrix inverse operation and found
some strange things. See the following results (the package Matrix is most
updated).
1. load the function Toeplitz and library(Matrix)
"Toeplitz"
function(x)
{
matrix(x[1 + abs(outer(seq(along = x), seq(along = x), FUN =
"-"))],
byrow = T, ncol = length(x))
}
2. On a PIII 1.125GHz (256kb cache)desktop PC with win2k, the results
are> system.time(solve(Toeplitz(.5^seq(0,999))))
[1] 38.25 0.22 42.01 NA NA> system.time(solve.Matrix(Toeplitz(.5^seq(0,999))))
[1] 28.01 0.17 30.94 NA NA
On a PIII-M 1.2GHz (512kb cache) laptop PC with winxp pro, the results
are> system.time(solve(Toeplitz(.5^seq(0,999))))
[1] 33.94 0.25 34.01 NA NA> system.time(solve.Matrix(Toeplitz(.5^seq(0,999))))
[1] 9.40 0.13 9.54 NA NA
On a PIII Xeon 1GHz (256kb cache) RedHat 7.2 Linux, the results
are> system.time(solve(Toeplitz(.5^seq(0,999))))
[1] 53.79 0.48 00.01 NA NA> system.time(solve.Matrix(Toeplitz(.5^seq(0,999))))
[1] 30.29 0.21 30.59 NA NA
This was the first time that windows R outperformed linux R even PIII Xeon 1GHz
is slightly faster than PIII 1.125GHz. Notice the abnormal result 9.4s from
PIII-M 1.2GHz (cannot be three times faster). After using native fortran lapack
(call to dpotrf and dpotri returns 9s for PIII Xeon), I realized that R may
treat the matrix as a general matrix rather than symmetric except the abnormal
case. To confirm this, I modified the R_LapackPP.cc in Matrix package and it
turned out that
is.MMatrix() failed to return true.
Hence
if (checkClass(classes, "Hermitian"))
return new LaSymmMatDouble(x);
is not executed in static LaMatDouble* asLaMatrix(SEXP x). Instead
if (isMatrix(x)) {
return new LaGenMatDouble(x);
}
is used.
After replaced
return new LaGenMatDouble(x);
with
return new LaSymmMatDouble(x);
the result was (PIII Xeon 1GHz RedHat 7.2 Linux))> system.time(solve.Matrix(Toeplitz(.5^seq(0,999))))
[1] 9.76 0.22 9.97 0.00 0.00
Clearly there are some bugs in solve() and solve.Matrix(). Notice that all
have> is.Hermitian(Toeplitz(.5^seq(0,999)))
[1] TRUE
Thanks
Hao Yu
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-devel mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !) To:
r-devel-request@stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._