This subroutine uses Householder transformations with column pivoting (optional) to compute a QR factorization of the m by n matrix A. That is, qrfac determines an orthogonal matrix Q, a permutation matrix P, and an upper trapezoidal matrix R with diagonal elements of nonincreasing magnitude, such that AP = QR. The Householder transformation for column k, k = 1,2,...,min(m,n), is of the form