From: <rvaradhan_at_jhmi.edu>

Date: Wed, 17 Jun 2009 23:45:11 +0200 (CEST)

https://stat.ethz.ch/mailman/listinfo/r-devel Received on Thu 18 Jun 2009 - 07:02:34 GMT

Date: Wed, 17 Jun 2009 23:45:11 +0200 (CEST)

Full_Name: Ravi Varadhan

Version: 2.8.1

OS: Windows

Submission from: (NULL) (162.129.251.19)

Inverting a matrix with solve(), but using LAPACK=TRUE, gives erroneous results:

hilbert <- function(n) { i <- 1:n; 1 / outer(i - 1, i, "+") } h5 <- hilbert(5) hinv1 <- solve(qr(h5)) hinv2 <- solve(qr(h5, LAPACK=TRUE)) all.equal(hinv1, hinv2) # They are not equal

Here is a function that I wrote to correct this problem:

solve.lapack <- function(A, LAPACK=TRUE, tol=1.e-07) { # A function to invert a matrix using "LAPACK" or "LINPACK" if (nrow(A) != ncol(A)) stop("Matrix muxt be square") qrA <- qr(A, LAPACK=LAPACK, tol=tol) if (LAPACK) { apply(diag(1, ncol(A)), 2, function(x) solve(qrA, x)) } else solve(qrA) } hinv3 <- solve.lapack(h5) all.equal(hinv1, hinv3) # Now, they are equal ______________________________________________R-devel_at_r-project.org mailing list

https://stat.ethz.ch/mailman/listinfo/r-devel Received on Thu 18 Jun 2009 - 07:02:34 GMT

Archive maintained by Robert King, hosted by
the discipline of
statistics at the
University of Newcastle,
Australia.

Archive generated by hypermail 2.2.0, at Thu 18 Jun 2009 - 13:36:25 GMT.

*
Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-devel.
Please read the posting
guide before posting to the list.
*