[Rd] Matrix identification bug (PR#1361)

[hyu@stats.uwo.ca]

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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._