Re: [Rd] NaN and linear algebra

From: Kjell Konis <konis_at_stats.ox.ac.uk>
Date: Thu 24 Mar 2005 - 20:03:45 GMT

> If, as reported, the behaviour is fundamentally similar on MacOS X, 
> Linux and IRIX then it would appear that this is a feature(bug?) of 
> the LAPACK reference implementation and that whoever wrote the Windows 
> version saw fit to improve on it.  I don't think LAPACK ever pretended 
> to properly handle NaNs.

Although both Linux and OSX produce an error the behavior of LAPACK is not similar.

(It would be nice if someone using Windows could try the following test too)

(1) Compile the following C code

void testdgesv(long* N, long* NRHS, double* A, long* LDA, long* IPIVOT, double* B, long* LDB, long* INFO)
{

   DGESV(N, NRHS, A, LDA, IPIVOT, B, LDB, INFO); }

(you will have to use the right combination of case and _s)

(2) Run these commands in R:

dyn.load("test.so")
N <- as.integer(3)
NRHS <- as.integer(3)
A <- matrix(as.double(NaN), 3, 3)
LDA <- as.integer(3)
IPIVOT <- integer(3)
B <- diag(3)
storage.mode(B) <- "double"
LDB <- as.integer(3)
INFO <- as.integer(0)

.C("testdgesv", N = N, NRHS = NRHS, A = A, LDA = LDA, IPIVOT = IPIVOT, B = B, LDB = LDB, INFO = INFO, NAOK = TRUE) $N
[1] 3

$NRHS
[1] 3

$A

      [,1] [,2] [,3]
[1,] NaN NaN NaN
[2,] NaN NaN NaN
[3,] NaN NaN NaN

$LDA
[1] 3

$IPIVOT
[1] 1 2 3

$B

      [,1] [,2] [,3]
[1,] 1 0 0
[2,] 0 1 0
[3,] 0 0 1

$LDB
[1] 3

$INFO
[1] 3

//////////////////

On my Mac mini this exits with no errors or warnings (other than INFO = 3 anyway). However,

solve(A)

generates the error:

Error in solve.default(A) : Lapack routine dgesv: system is exactly singular

This behavior seems to suggest that error is generated somewhere in R, not in LAPACK.

Okay, now for the weird bit.

When I run this in Linux I get:

 > .C("testdgesv", N = N, NRHS = NRHS, A = A, LDA = LDA, IPIVOT = IPIVOT, B = B, LDB = LDB, INFO = INFO, NAOK = TRUE) $N
[1] 3

$NRHS
[1] 3

$A

      [,1] [,2] [,3]
[1,] NaN NaN NaN
[2,] NaN NaN NaN
[3,] NaN NaN NaN

$LDA
[1] 3

$IPIVOT
[1] 3 2 3

$B

      [,1] [,2] [,3]
[1,] NaN NaN NaN
[2,] NaN NaN NaN
[3,] NaN NaN NaN

$LDB
[1] 3

$INFO
[1] 0

where INFO == 0 indicates "successful exit" but

 > solve(A)
Error in solve.default(A) : system is computationally singular: reciprocal condition number = 0

In case it is important:

platform powerpc-apple-darwin7.8.0

arch     powerpc
os       darwin7.8.0

system powerpc, darwin7.8.0
status
major 2
minor 0.1
year 2004
month 11
day 15
language R

and

platform i686-pc-linux-gnu

arch     i686
os       linux-gnu

system i686, linux-gnu
status
major 2
minor 0.1
year 2004
month 11
day 15
language R

prompt> more Makevars
PKG_LIBS=${LAPACK_LIBS} ${BLAS_LIBS} ${FLIBS}



R-devel@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel Received on Fri Mar 25 07:11:32 2005

This archive was generated by hypermail 2.1.8 : Mon 24 Oct 2005 - 22:26:32 GMT