[Rd] Inverting a square matrix using solve() with LAPACK=TRUE (PR#13762)

rvaradhan at jhmi.edu rvaradhan at jhmi.edu
Wed Jun 17 23:45:11 CEST 2009

Full_Name: Ravi Varadhan
Version: 2.8.1
OS: Windows
Submission from: (NULL) (

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

Here is an example:

     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

More information about the R-devel mailing list