dyn.load ("matrix.so") gsRC <- function (x, eps = 1e-10) { n <- nrow (x) m <- ncol (x) h <- .C( "gsC", x = as.double(x), r = as.double (matrix (0, m, m)), n = as.integer (n), m = as.integer (m), rank = as.integer (0), pivot = as.integer (1:m), eps = as.double (eps) ) rank = h\$rank return (list ( q = matrix (h\$x, n, m)[, 1:rank, drop = FALSE], r = matrix (h\$r, m, m)[1:rank, , drop = FALSE], rank = rank, pivot = h\$pivot )) } lsRC <- function (x, y, eps = 1e-10) { x <- as.matrix (x) y <- as.matrix (y) n <- nrow (x) m <- ncol (x) h <- gsRC (x, eps) l <- h\$rank p <- order (h\$pivot) k <- 1:l q <- h\$q a <- h\$r[, k, drop = FALSE] v <- h\$r[, -k, drop = FALSE] u <- crossprod (q, y) b <- solve (a, u) res <- drop (y - q %*% u) s <- sum (res ^ 2) b <- rbind(b, matrix (0, m - l, ncol(y)))[p, , drop = FALSE] if (l == m) { e <- matrix(0, m, 1) } else { e <- rbind (-solve(a, v), diag(m - l))[p, , drop = FALSE] } return (list ( solution = b, residuals = res, minssq = s, nullspace = e, rank = l, pivot = p )) } nullRC <- function (x, eps = 1e-10) { h <- gsRC (x, eps = eps) rank <- h\$rank r <- h\$r m <- ncol (x) t <- r[, 1:rank, drop = FALSE] s <- r[, -(1:rank), drop = FALSE] if (rank == m) return (matrix(0, m, 1)) else { nullspace <- rbind (-solve(t, s), diag (m - rank))[order(h\$pivot), , drop = FALSE] return (gsRC (nullspace)\$q) } } ginvRC <- function (x, eps = 1e-10) { h <- gsRC (x, eps) p <- order(h\$pivot) q <- h\$q s <- h\$r z <- crossprod (s, (solve (tcrossprod(s), t(q)))) return (z[p, , drop = FALSE]) }