From 66761e1f3787d6a3ea297260136901e969ac9e14 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jiahao=20Chen=20=28=E9=99=88=E5=AE=B6=E8=B1=AA=29?=
 <jiahao@mit.edu>
Date: Wed, 23 Oct 2013 19:20:22 -0400
Subject: [PATCH 1/8] Refactors QR code

- Deprecates qrp, qrpfact and qrpfact!
- Adds trailing pivoted::Bool option to qr, qrfact and qrpfact!
- Reorganizes QR factorizations section so that methods common to QR,
  QRPackedQ, QRPivoted, QRPivotedQ are grouped together.
- Remove improperly cleaned up blob from merge
- Replace methods defining default parameters with methods with default
  parameter values
- Rename QRPivoted.jpvt to QRPivoted.pivots
---
 base/deprecated.jl           |   4 +
 base/linalg/dense.jl         |   4 +-
 base/linalg/factorization.jl | 233 ++++++++++++++++-------------------
 3 files changed, 109 insertions(+), 132 deletions(-)

diff --git a/base/deprecated.jl b/base/deprecated.jl
index fce51d43c4e88..98ddbbc4eba4f 100644
--- a/base/deprecated.jl
+++ b/base/deprecated.jl
@@ -168,6 +168,10 @@ export PipeString
 #   copy!(dest::AbstractArray, doffs::Integer, src::Integer)
 # in abstractarray.jl
 @deprecate copy!(dest::AbstractArray, src, doffs::Integer)  copy!(dest, doffs, src)
+@deprecate qrpfact!(A)         qrfact!(A, true)
+@deprecate qrpfact(A)          qrfact(A, true)
+@deprecate qrp(A, thin)        qr(A, thin, true)
+@deprecate qrp(A)              qr(A, true)
 
 deprecated_ls() = run(`ls -l`)
 deprecated_ls(args::Cmd) = run(`ls -l $args`)
diff --git a/base/linalg/dense.jl b/base/linalg/dense.jl
index aaa2a29ba9b41..6034132055d04 100644
--- a/base/linalg/dense.jl
+++ b/base/linalg/dense.jl
@@ -496,7 +496,7 @@ function factorize!{T}(A::Matrix{T})
         end
         return lufact!(A)
     end
-    return qrpfact!(A)
+    return qrfact!(A, true)
 end
 
 factorize(A::AbstractMatrix) = factorize!(copy(A))
@@ -517,7 +517,7 @@ function (\){T<:BlasFloat}(A::StridedMatrix{T}, B::StridedVecOrMat{T})
         if istriu(A) return \(Triangular(A, :U),B) end
         return \(lufact(A),B)
     end
-    return qrpfact(A)\B
+    return qrfact(A, true)\B
 end
 
 ## Moore-Penrose inverse
diff --git a/base/linalg/factorization.jl b/base/linalg/factorization.jl
index f6c56022f40c5..269434b70430c 100644
--- a/base/linalg/factorization.jl
+++ b/base/linalg/factorization.jl
@@ -244,69 +244,112 @@ end
 
 cond(A::LU, p) = 1.0/LinAlg.LAPACK.gecon!(p == 1 ? '1' : 'I', A.factors, norm(A[:L][A[:p],:]*A[:U], p))
 
-## QR decomposition without column pivots. By the faster geqrt3
+####################
+# QR factorization #
+####################
+
 type QR{S<:BlasFloat} <: Factorization{S}
     vs::Matrix{S}                     # the elements on and above the diagonal contain the N-by-N upper triangular matrix R; the elements below the diagonal are the columns of V
     T::Matrix{S}                      # upper triangular factor of the block reflector.
 end
+
+# QR decomposition without column pivots. By the faster geqrt3
 QR{T<:BlasFloat}(A::StridedMatrix{T}, nb::Integer = min(minimum(size(A)), 36)) = QR(LAPACK.geqrt!(A, nb)...)
 
-qrfact!{T<:BlasFloat}(A::StridedMatrix{T}, args::Integer...) = QR(A, args...)
-qrfact!(A::StridedMatrix, args::Integer...) = qrfact!(float(A), args...)
-qrfact{T<:BlasFloat}(A::StridedMatrix{T}, args::Integer...) = qrfact!(copy(A), args...)
-qrfact(A::StridedMatrix, args::Integer...) = qrfact!(float(A), args...)
-qrfact(x::Integer) = qrfact(float(x))
-qrfact(x::Number) = QR(fill(one(x), 1, 1), fill(x, 1, 1))
+type QRPackedQ{S} <: AbstractMatrix{S} 
+    vs::Matrix{S}
+    T::Matrix{S}
+end
+QRPackedQ(A::QR) = QRPackedQ(A.vs, A.T)
 
-function qr(A::Union(Number, AbstractMatrix), thin::Bool)
-    F = qrfact(A)
-    return (full(F[:Q], thin), F[:R])
+type QRPivoted{T} <: Factorization{T}
+    hh::Matrix{T}
+    tau::Vector{T}
+    pivots::Vector{BlasInt}
 end
-qr(A::Union(Number, AbstractMatrix)) = qr(A, true)
 
-size(A::QR, args::Integer...) = size(A.vs, args...)
+type QRPivotedQ{T} <: AbstractMatrix{T}
+    hh::Matrix{T}                      # Householder transformations and R
+    tau::Vector{T}                     # Scalar factors of transformations
+end
+QRPivotedQ(A::QRPivoted) = QRPivotedQ(A.hh, A.tau)
 
-function getindex(A::QR, d::Symbol)
-    if d == :R; return triu(A.vs[1:minimum(size(A)),:]); end;
-    if d == :Q; return QRPackedQ(A); end
-    error("No such type field")
+function qr(A::Union(Number, AbstractMatrix), thin::Bool=true, pivoted::Bool=false)
+  F = qrfact(A, pivoted)
+  return pivoted ? (full(F[:Q], thin), F[:R], F[:p]) : (full(F[:Q], thin), F[:R])
 end
 
-type QRPackedQ{S} <: AbstractMatrix{S} 
-    vs::Matrix{S}                      
-    T::Matrix{S}                       
+qrfact!{T<:BlasFloat}(A::StridedMatrix{T}, pivoted::Bool=false, args::Integer...) = pivoted ? QRPivoted{T}(LAPACK.geqp3!(A)...) : QR(A, args...) 
+qrfact!(A::StridedMatrix, pivoted::Bool=false) = qrfact!(float(A), pivoted)
+qrfact{T<:BlasFloat}(A::StridedMatrix{T}, pivoted::Bool=false, args::Integer...) = qrfact!(copy(A), pivoted, args...)
+qrfact(A::StridedMatrix, pivoted::Bool=false, args::Integer...) = qrfact!(float(A), pivoted, args...)
+qrfact(x::Integer, pivoted::Bool=false) = qrfact(float(x), pivoted)
+qrfact(x::Number, pivoted::Bool=false) = pivoted ? error("Not implemented") : QR(fill(one(x), 1, 1), fill(x, 1, 1))
+size(A::Union(QR,        QRPackedQ ), args::Integer...) = size(A.vs, args...)
+size(A::Union(QRPivoted, QRPivotedQ), args::Integer...) = size(A.hh, args...)
+
+function getindex(A::QR, d::Symbol)
+    if d == :R return triu(A.vs[1:minimum(size(A)),:])
+    elseif d == :Q return QRPackedQ(A)
+    else error("No such type field") end
 end
-QRPackedQ(A::QR) = QRPackedQ(A.vs, A.T)
 
-size(A::QRPackedQ, args::Integer...) = size(A.vs, args...)
+function getindex{T<:BlasFloat}(A::QRPivoted{T}, d::Symbol)
+    if d == :R return triu(A.hh[1:minimum(size(A)),:])
+    elseif d == :Q return QRPivotedQ(A)
+    elseif d == :p return A.pivots
+    elseif d == :P #The permutation matrix for the pivot
+        p = A[:p]
+        n = length(p)
+        P = zeros(T, n, n)
+        for i in 1:n
+            P[p[i],i] = one(T)
+        end
+        return P
+    end
+    error("No such type field")
+end
 
-function full{T<:BlasFloat}(A::QRPackedQ{T}, thin::Bool)
-    if thin return A * eye(T, size(A.T, 2)) end
-    return A * eye(T, size(A, 1))
+full{T<:BlasFloat}(A::QRPackedQ{T}, thin::Bool) = A * eye(T, thin ? size(A.T, 2) : size(A, 1))
+function full{T<:BlasFloat}(A::QRPivotedQ{T}, thin::Bool=true)
+    m, n = size(A.hh)
+    return LAPACK.orgqr!(thin ? copy(A.hh) : [A.hh zeros(T, m, max(0, m-n))], A.tau)
 end
-full(A::QRPackedQ) = full(A, true)
 
-print_matrix(io::IO, A::QRPackedQ, rows::Integer, cols::Integer) = print_matrix(io, full(A), rows, cols)
+print_matrix(io::IO, A::Union(QRPackedQ, QRPivotedQ), rows::Integer, cols::Integer) = print_matrix(io, full(A), rows, cols)
 
 ## Multiplication by Q from the QR decomposition
 function *{T<:BlasFloat}(A::QRPackedQ{T}, B::StridedVecOrMat{T})
-    m = size(B, 1)
-    n = size(B, 2)
+    m, n = size(B)
     if m == size(A.vs, 1)
         Bc = copy(B)
     elseif m == size(A.vs, 2)
-        Bc = [B; zeros(T, size(A.vs, 1) - m, n)]
+        Bc = [B; zeros(T, size(A.vs,1)-m, n)]
     else
         throw(DimensionMismatch(""))
     end
     LAPACK.gemqrt!('L', 'N', A.vs, A.T, Bc)
 end
-Ac_mul_B{T<:BlasReal}(A::QRPackedQ{T}, B::StridedVecOrMat{T}) = LAPACK.gemqrt!('L','T',A.vs,A.T,copy(B))
-Ac_mul_B{T<:BlasComplex}(A::QRPackedQ{T}, B::StridedVecOrMat{T}) = LAPACK.gemqrt!('L','C',A.vs,A.T,copy(B))
 *{T<:BlasFloat}(A::StridedVecOrMat{T}, B::QRPackedQ{T}) = LAPACK.gemqrt!('R', 'N', B.vs, B.T, copy(A))
+function *{T<:BlasFloat}(A::QRPivotedQ{T}, B::StridedVecOrMat{T})
+    m, n = size(B)
+    if m == size(A.hh, 1)
+        Bc = copy(B)
+    elseif m == size(A.hh, 2)
+        Bc = [B; zeros(T, size(A.hh,1)-m, n)]
+    else
+        throw(DimensionMismatch(""))
+    end
+    LAPACK.ormqr!('L', 'N', A.hh, A.tau, Bc)
+end
+*(A::StridedVecOrMat, B::QRPivotedQ) = LAPACK.ormqr!('R', 'N', B.hh, B.tau, copy(A))
+
+Ac_mul_B{T<:BlasReal   }(A::QRPackedQ{T},  B::StridedVecOrMat) = LAPACK.gemqrt!('L','T',A.vs,A.T,  copy(B))
+Ac_mul_B{T<:BlasComplex}(A::QRPackedQ{T},  B::StridedVecOrMat) = LAPACK.gemqrt!('L','C',A.vs,A.T,  copy(B))
+Ac_mul_B{T<:BlasReal   }(A::QRPivotedQ{T}, B::StridedVecOrMat) = LAPACK.ormqr! ('L','T',A.hh,A.tau,copy(B))
+Ac_mul_B{T<:BlasComplex}(A::QRPivotedQ{T}, B::StridedVecOrMat) = LAPACK.ormqr! ('L','C',A.hh,A.tau,copy(B))
 function A_mul_Bc{T<:BlasFloat}(A::StridedVecOrMat{T}, B::QRPackedQ{T})
-    m = size(A, 1)
-    n = size(A, 2)
+    m, n = size(A)
     if n == size(B.vs, 1)
         Ac = copy(A)
     elseif n == size(B.vs, 2)
@@ -316,127 +359,57 @@ function A_mul_Bc{T<:BlasFloat}(A::StridedVecOrMat{T}, B::QRPackedQ{T})
     end
     LAPACK.gemqrt!('R', iseltype(B.vs,Complex) ? 'C' : 'T', B.vs, B.T, Ac)
 end
-## Least squares solution.  Should be more careful about cases with m < n
-(\)(A::QR, B::StridedVector) = Triangular(A[:R], :U)\(A[:Q]'B)[1:size(A, 2)]
-(\)(A::QR, B::StridedMatrix) = Triangular(A[:R], :U)\(A[:Q]'B)[1:size(A, 2),:]
-
-type QRPivoted{T} <: Factorization{T}
-    hh::Matrix{T}
-    tau::Vector{T}
-    jpvt::Vector{BlasInt}
-end
-
-qrpfact!{T<:BlasFloat}(A::StridedMatrix{T}) = QRPivoted{T}(LAPACK.geqp3!(A)...)
-qrpfact!(A::StridedMatrix) = qrpfact!(float(A))
-qrpfact{T<:BlasFloat}(A::StridedMatrix{T}) = qrpfact!(copy(A))
-qrpfact(A::StridedMatrix) = qrpfact!(float(A))
-
-function qrp(A::AbstractMatrix, thin::Bool)
-    F = qrpfact(A)
-    return full(F[:Q], thin), F[:R], F[:p]
-end
-qrp(A::AbstractMatrix) = qrp(A, false)
-
-size(A::QRPivoted, args::Integer...) = size(A.hh, args...)
-
-function getindex{T<:BlasFloat}(A::QRPivoted{T}, d::Symbol)
-    if d == :R; return triu(A.hh[1:minimum(size(A)),:]); end;
-    if d == :Q; return QRPivotedQ(A); end
-    if d == :p; return A.jpvt; end
-    if d == :P
-        p = A[:p]
-        n = length(p)
-        P = zeros(T, n, n)
-        for i in 1:n
-            P[p[i],i] = one(T)
-        end
-        return P
+function A_mul_Bc{T<:BlasFloat}(A::StridedVecOrMat{T}, B::QRPivotedQ{T})
+    m, n = size(A)
+    if n == size(B.hh, 1)
+        Ac = copy(A)
+    elseif n == size(B.hh, 2)
+        Ac = [B zeros(T, m, size(B.hh, 1) - n)]
+    else
+        throw(DimensionMismatch(""))
     end
-    error("No such type field")
+    LAPACK.ormqr!('R', iseltype(B.hh,Complex) ? 'C' : 'T', B.hh, B.tau, Ac)
 end
 
+## Least squares solution.
+#XXX Should be more careful about cases with m < n
+(\)(A::QR, B::StridedVector) = Triangular(A[:R], :U)\(A[:Q]'B)[1:size(A, 2)]
+(\)(A::QR, B::StridedMatrix) = Triangular(A[:R], :U)\(A[:Q]'B)[1:size(A, 2),:]
 # Julia implementation similarly to xgelsy
 function (\){T<:BlasFloat}(A::QRPivoted{T}, B::StridedMatrix{T}, rcond::Real)
-    nr = minimum(size(A.hh))
-    nrhs = size(B, 2)
+    nr, nrhs = minimum(size(A.hh)), size(B, 2)
     if nr == 0 return zeros(0, nrhs), 0 end
     ar = abs(A.hh[1])
-    if ar == 0 return zeros(nr, nrhs), 0 end
+    if ar < eps(typeof(ar)) return zeros(nr, nrhs), 0 end
     rnk = 1
-    xmin = ones(T, nr)
-    xmax = ones(T, nr)
-    tmin = ar
-    tmax = ar
+    xmin, xmax = ones(T, nr), ones(T, nr)
+    tmin, tmax = ar, ar
     while rnk < nr
-        tmin, smin, cmin = LAPACK.laic1!(2, sub(xmin, 1:rnk), tmin, sub(A.hh, 1:rnk, rnk + 1), A.hh[rnk + 1, rnk + 1])
-        tmax, smax, cmax = LAPACK.laic1!(1, sub(xmax, 1:rnk), tmax, sub(A.hh, 1:rnk, rnk + 1), A.hh[rnk + 1, rnk + 1])
+        Ahh = sub(A.hh, 1:rnk, rnk+1), A.hh[rnk+1, rnk+1]
+        tmin, smin, cmin = LAPACK.laic1!(2, sub(xmin, 1:rnk), tmin, Ahh...)
+        tmax, smax, cmax = LAPACK.laic1!(1, sub(xmax, 1:rnk), tmax, Ahh...)
         if tmax*rcond > tmin break end
-        xmin[1:rnk + 1] = [smin*sub(xmin, 1:rnk), cmin]
-        xmax[1:rnk + 1] = [smax*sub(xmin, 1:rnk), cmax]
+        xmin[1:rnk+1] = [smin*sub(xmin, 1:rnk), cmin]
+        xmax[1:rnk+1] = [smax*sub(xmin, 1:rnk), cmax] #XXX should this be xmax?
         rnk += 1
         # if cond(r[1:rnk, 1:rnk])*rcond < 1 break end
     end
     C, tau = LAPACK.tzrzf!(A.hh[1:rnk,:])
-    X = [Triangular(C[1:rnk,1:rnk],:U)\(A[:Q]'B)[1:rnk,:]; zeros(T, size(A.hh, 2) - rnk, nrhs)]
+    X = [Triangular(C[1:rnk,1:rnk],:U)\(A[:Q]'B)[1:rnk,:]; zeros(T, size(A.hh,2)-rnk, nrhs)]
     LAPACK.ormrz!('L', iseltype(B, Complex) ? 'C' : 'T', C, tau, X)
     return X[invperm(A[:p]),:], rnk
 end
 (\)(A::QRPivoted, B::StridedMatrix) = (\)(A, B, sqrt(eps(typeof(real(B[1])))))[1]
 (\)(A::QRPivoted, B::StridedVector) = (\)(A, reshape(B, length(B), 1))[:]
 
-type QRPivotedQ{T} <: AbstractMatrix{T}
-    hh::Matrix{T}                       # Householder transformations and R
-    tau::Vector{T}                      # Scalar factors of transformations
-end
-QRPivotedQ(A::QRPivoted) = QRPivotedQ(A.hh, A.tau)
-
-size(A::QRPivotedQ, args...) = size(A.hh, args...)
-
-function full{T<:BlasFloat}(A::QRPivotedQ{T}, thin::Bool)
-    m, n = size(A.hh)
-    if !thin
-        B = [A.hh zeros(T, m, max(0, m - n))]
-        return LAPACK.orgqr!(B, A.tau)
-    end
-    return LAPACK.orgqr!(copy(A.hh), A.tau)
-end
-
-full(A::QRPivotedQ) = full(A, true)
-print_matrix(io::IO, A::QRPivotedQ, rows::Integer, cols::Integer) = print_matrix(io, full(A), rows, cols)
-
-## Multiplication by Q from the Pivoted QR decomposition
-function *{T<:BlasFloat}(A::QRPivotedQ{T}, B::StridedVecOrMat{T})
-    m = size(B, 1)
-    n = size(B, 2)
-    if m == size(A.hh, 1)
-        Bc = copy(B)
-    elseif m == size(A.hh, 2)
-        Bc = [B; zeros(T, size(A.hh, 1) - m, n)]
-    else
-        throw(DimensionMismatch(""))
-    end
-    LAPACK.ormqr!('L', 'N', A.hh, A.tau, Bc)
-end
-Ac_mul_B{T<:BlasReal}(A::QRPivotedQ{T}, B::StridedVecOrMat{T}) = LAPACK.ormqr!('L','T',A.hh,A.tau,copy(B))
-Ac_mul_B{T<:BlasComplex}(A::QRPivotedQ{T}, B::StridedVecOrMat{T}) = LAPACK.ormqr!('L','C',A.hh,A.tau,copy(B))
-*(A::StridedVecOrMat, B::QRPivotedQ) = LAPACK.ormqr!('R', 'N', B.hh, B.tau, copy(A))
-function A_mul_Bc{T<:BlasFloat}(A::StridedVecOrMat{T}, B::QRPivotedQ{T})
-    m = size(A, 1)
-    n = size(A, 2)
-    if n == size(B.hh, 1)
-        Ac = copy(A)
-    elseif n == size(B.hh, 2)
-        Ac = [B zeros(T, m, size(B.hh, 1) - n)]
-    else
-        throw(DimensionMismatch(""))
-    end
-    LAPACK.ormqr!('R', iseltype(B.hh,Complex) ? 'C' : 'T', B.hh, B.tau, Ac)
-end
-
 ##TODO:  Add methods for rank(A::QRP{T}) and adjust the (\) method accordingly
 ##       Add rcond methods for Cholesky, LU, QR and QRP types
 ## Lower priority: Add LQ, QL and RQ factorizations
 
+############################
+# Hessenberg factorization #
+############################
+
 # FIXME! Should add balancing option through xgebal
 type Hessenberg{T} <: Factorization{T}
     hh::Matrix{T}

From 6348a3a3963c6585edef0025303445a57b7fe1ab Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jiahao=20Chen=20=28=E9=99=88=E5=AE=B6=E8=B1=AA=29?=
 <jiahao@mit.edu>
Date: Thu, 24 Oct 2013 17:38:36 -0400
Subject: [PATCH 2/8] Deprecates cholpfact, cholpfact!

---
 base/deprecated.jl           |   4 +
 base/linalg/factorization.jl | 197 +++++++++++++++--------------------
 2 files changed, 90 insertions(+), 111 deletions(-)

diff --git a/base/deprecated.jl b/base/deprecated.jl
index 98ddbbc4eba4f..25b5a7facdc58 100644
--- a/base/deprecated.jl
+++ b/base/deprecated.jl
@@ -172,6 +172,10 @@ export PipeString
 @deprecate qrpfact(A)          qrfact(A, true)
 @deprecate qrp(A, thin)        qr(A, thin, true)
 @deprecate qrp(A)              qr(A, true)
+@deprecate cholpfact!(A)            cholfact!(A, :U, true)
+@deprecate cholpfact!(A,tol)        cholfact!(A, :U, true, tol)
+@deprecate cholpfact!(A,uplo,tol)   cholfact!(A, uplo, true, tol)
+@deprecate cholpfact(A)             cholfact(A, :U, true)
 
 deprecated_ls() = run(`ls -l`)
 deprecated_ls(args::Cmd) = run(`ls -l $args`)
diff --git a/base/linalg/factorization.jl b/base/linalg/factorization.jl
index 269434b70430c..2a95dda5548c0 100644
--- a/base/linalg/factorization.jl
+++ b/base/linalg/factorization.jl
@@ -4,61 +4,15 @@ abstract Factorization{T}
 
 (\)(F::Factorization, b::Union(AbstractVector, AbstractMatrix)) = A_ldiv_B!(F, copy(b))
 
+##########################
+# Cholesky factorization #
+##########################
+
 type Cholesky{T<:BlasFloat} <: Factorization{T}
     UL::Matrix{T}
     uplo::Char
 end
 
-function cholfact!{T<:BlasFloat}(A::StridedMatrix{T}, uplo::Symbol)
-    uplochar = string(uplo)[1]
-    C, info = LAPACK.potrf!(uplochar, A)
-    if info > 0 throw(PosDefException(info)) end
-    Cholesky(C, uplochar)
-end
-cholfact!(A::StridedMatrix, args...) = cholfact!(float(A), args...)
-cholfact!{T<:BlasFloat}(A::StridedMatrix{T}) = cholfact!(A, :U)
-cholfact{T<:BlasFloat}(A::StridedMatrix{T}, args...) = cholfact!(copy(A), args...)
-cholfact(A::StridedMatrix, args...) = cholfact!(float(A), args...)
-cholfact(x::Number) = imag(x) == 0 && real(x) > 0 ? Cholesky(fill(sqrt(x), 1, 1), :U) : throw(PosDefException(1))
-
-chol(A::Union(Number, AbstractMatrix), uplo::Symbol) = cholfact(A, uplo)[uplo]
-chol(A::Union(Number, AbstractMatrix)) = cholfact(A, :U)[:U]
-
-size(C::Cholesky) = size(C.UL)
-size(C::Cholesky,d::Integer) = size(C.UL,d)
-
-function getindex(C::Cholesky, d::Symbol)
-    C.uplo == 'U' ? triu!(C.UL) : tril!(C.UL)
-    if d == :U || d == :L
-        return symbol(C.uplo) == d ? C.UL : C.UL'
-    elseif d == :UL
-        return Triangular(C.UL, C.uplo)
-    end
-    error("No such type field")
-end
-
-A_ldiv_B!{T<:BlasFloat}(C::Cholesky{T}, B::StridedVecOrMat{T}) =
-    LAPACK.potrs!(C.uplo, C.UL, B)
-
-function det{T}(C::Cholesky{T})
-    dd = one(T)
-    for i in 1:size(C.UL,1) dd *= abs2(C.UL[i,i]) end
-    dd
-end
-
-function logdet{T}(C::Cholesky{T})
-    dd = zero(T)
-    for i in 1:size(C.UL,1) dd += log(C.UL[i,i]) end
-    dd + dd # instead of 2.0dd which can change the type
-end
-
-function inv(C::Cholesky)
-    Ci, info = LAPACK.potri!(C.uplo, copy(C.UL))
-    if info != 0; throw(SingularException(info)); end 
-    symmetrize_conj!(Ci, C.uplo)
-end
-
-## Pivoted Cholesky
 type CholeskyPivoted{T<:BlasFloat} <: Factorization{T}
     UL::Matrix{T}
     uplo::Char
@@ -67,44 +21,70 @@ type CholeskyPivoted{T<:BlasFloat} <: Factorization{T}
     tol::Real
     info::BlasInt
 end
+
 function CholeskyPivoted{T<:BlasFloat}(A::StridedMatrix{T}, uplo::Char, tol::Real)
     A, piv, rank, info = LAPACK.pstrf!(uplo, A, tol)
-    CholeskyPivoted{T}(uplo == 'U' ? triu!(A) : tril!(A), uplo, piv, rank, tol, info)
+    CholeskyPivoted{T}(uplo=='U' ? triu!(A) : tril!(A), uplo, piv, rank, tol, info)
+end
+
+function cholfact!{T<:BlasFloat}(A::StridedMatrix{T}, uplo::Symbol=:U, pivoted::Bool=false,
+      tol::Real=-1.0)
+  if pivoted
+    CholeskyPivoted(A, string(uplo)[1], tol)
+  else
+    uplochar = string(uplo)[1]
+    C, info = LAPACK.potrf!(uplochar, A)
+    info==0 ? Cholesky(C, uplochar) : throw(PosDefException(info))
+  end
 end
-cholpfact!(A::StridedMatrix, args...) = cholpfact!(float(A), args...)
-cholpfact!{T<:BlasFloat}(A::StridedMatrix{T}, uplo::Symbol, tol::Real) = CholeskyPivoted(A, string(uplo)[1], tol)
-cholpfact!{T<:BlasFloat}(A::StridedMatrix{T}, tol::Real) = cholpfact!(A, :U, tol)
-cholpfact!{T<:BlasFloat}(A::StridedMatrix{T}) = cholpfact!(A, -1.)
-cholpfact{T<:BlasFloat}(A::StridedMatrix{T}, args...) = cholpfact!(copy(A), args...)
-cholpfact(A::StridedMatrix, args...) = cholpfact!(float(A), args...)
 
-size(C::CholeskyPivoted) = size(C.UL)
-size(C::CholeskyPivoted,d::Integer) = size(C.UL,d)
+cholfact!(A::StridedMatrix, uplo::Symbol=:U, pivoted::Bool=false, args...) = cholfact!(float(A), uplo, pivoted, args...)
+cholfact{T<:BlasFloat}(A::StridedMatrix{T}, uplo::Symbol=:U, pivoted::Bool=false, args...) = cholfact!(copy(A), uplo, pivoted, args...)
+cholfact(A::StridedMatrix, uplo::Symbol=:U, pivoted::Bool=false, args...) = cholfact!(float(A), uplo, pivoted, args...)
+cholfact(A::StridedMatrix, pivoted::Bool=false, args...) = cholfact(float(A), :U, pivoted, args...)
+cholfact(x::Number, pivoted::Bool=false) = imag(x) == 0 && real(x) > 0 ? Cholesky(fill(sqrt(x), 1, 1), :U) : throw(PosDefException(1))
 
+chol(A::Union(Number, AbstractMatrix), uplo::Symbol=:U) = cholfact(A, uplo)[uplo]
+
+function getindex(C::Cholesky, d::Symbol)
+    C.uplo == 'U' ? triu!(C.UL) : tril!(C.UL)
+    if d == :U || d == :L return symbol(C.uplo) == d ? C.UL : C.UL'
+    elseif d == :UL return Triangular(C.UL, C.uplo)
+    else error("No such type field")
+    end
+end
 getindex(C::CholeskyPivoted) = C.UL, C.piv
 function getindex{T<:BlasFloat}(C::CholeskyPivoted{T}, d::Symbol)
-    if d == :U || d == :L
-        return symbol(C.uplo) == d ? C.UL : C.UL'
+    if d == :U || d == :L return symbol(C.uplo) == d ? C.UL : C.UL'
+    elseif d == :p return C.piv 
+    elseif d == :P return perm_matrix(C.piv, T)
+    else error("No such type field")
     end
-    if d == :p return C.piv end
-    if d == :P
-        n = size(C, 1)
-        P = zeros(T, n, n)
-        for i in 1:n
-            P[C.piv[i],i] = one(T)
-        end
-        return P
+end
+
+#Generate permutation matrix from a permutation
+function perm_matrix{T<:Integer}(A::Vector{T}, output_type::Type=T)
+    n = length(A)
+    P = zeros(output_type, n, n)
+    for i=1:n
+        P[A[i], i] = one(output_type)
     end
-    error("No such type field")
+    P
 end
 
+size(C::Union(Cholesky,CholeskyPivoted)) = size(C.UL)
+size(C::Union(Cholesky,CholeskyPivoted),d::Integer) = size(C.UL,d)
+
+rank(C::CholeskyPivoted) = C.rank
+rank(C::Cholesky) = size(C.UL, 1) #full rank
+
+A_ldiv_B!{T<:BlasFloat}(C::Cholesky{T}, B::StridedVecOrMat{T}) = LAPACK.potrs!(C.uplo, C.UL, B)
 function A_ldiv_B!{T<:BlasFloat}(C::CholeskyPivoted{T}, B::StridedVector{T})
-    if C.rank < size(C.UL, 1); throw(RankDeficientException(C.info)); end
+    C.rank < size(C.UL, 1) ? throw(RankDeficientException(C.info)) :
     ipermute!(LAPACK.potrs!(C.uplo, C.UL, permute!(B, C.piv)), C.piv)
 end
-
 function A_ldiv_B!{T<:BlasFloat}(C::CholeskyPivoted{T}, B::StridedMatrix{T})
-    if C.rank < size(C.UL, 1); throw(RankDeficientException(C.info)); end
+    C.rank < size(C.UL, 1) ? throw(RankDeficientException(C.info)) : nothing
     n = size(C, 1)
     for i = 1:size(B, 2)
         permute!(sub(B, 1:n, i), C.piv)
@@ -113,19 +93,19 @@ function A_ldiv_B!{T<:BlasFloat}(C::CholeskyPivoted{T}, B::StridedMatrix{T})
     for i = 1:size(B, 2)
         ipermute!(sub(B, 1:n, i), C.piv)
     end
-    return B
+    B
 end
 
-rank(C::CholeskyPivoted) = C.rank
+det{T}(C::Cholesky{T}) = prod(abs2(diag(C.UL)))
+det{T}(C::CholeskyPivoted{T}) = C.rank < size(C.UL, 1) ? real(zero(T)) : prod(abs2(diag(C.UL)))
 
-function det{T}(C::CholeskyPivoted{T})
-    if C.rank < size(C.UL, 1) 
-        return real(zero(T))
-    else 
-        return prod(abs2(diag(C.UL)))
-    end
+logdet{T}(C::Cholesky{T}) = 2sum(log(diag(C.UL)))
+logdet{T}(C::CholeskyPivoted{T}) = C.rank < size(C.UL, 1) ? real(convert(T, -Inf)) : 2sum(log(diag(C.UL)))
+
+function inv(C::Cholesky)
+    Ci, info = LAPACK.potri!(C.uplo, copy(C.UL))
+    info == 0 ? symmetrize_conj!(Ci, C.uplo) : throw(SingularException(info))
 end
-    
 function inv(C::CholeskyPivoted)
     if C.rank < size(C.UL, 1) throw(RankDeficientException(C.info)) end
     Ci, info = LAPACK.potri!(C.uplo, copy(C.UL))
@@ -134,15 +114,18 @@ function inv(C::CholeskyPivoted)
     (symmetrize!(Ci, C.uplo))[ipiv, ipiv]
 end
 
-## LU
+####################
+# LU factorization #
+####################
+
 type LU{T<:BlasFloat} <: Factorization{T}
     factors::Matrix{T}
-    ipiv::Vector{BlasInt}
+    pivots::Vector{BlasInt}
     info::BlasInt
 end
 function LU{T<:BlasFloat}(A::StridedMatrix{T})
-    factors, ipiv, info = LAPACK.getrf!(A)
-    LU{T}(factors, ipiv, info)
+    factors, pivots, info = LAPACK.getrf!(A)
+    LU{T}(factors, pivots, info)
 end
 lufact!(A::StridedMatrix) = lufact!(float(A))
 lufact!{T<:BlasFloat}(A::StridedMatrix{T}) = LU(A)
@@ -160,18 +143,17 @@ size(A::LU,n) = size(A.factors,n)
 
 function getindex{T}(A::LU{T}, d::Symbol)
     m, n = size(A)
-    if d == :L; return tril(A.factors[1:m, 1:min(m,n)], -1) + eye(T, m, min(m,n)); end;
-    if d == :U; return triu(A.factors[1:min(m,n),1:n]); end;
-    if d == :p
+    if d == :L return tril(A.factors[1:m, 1:min(m,n)], -1) + eye(T, m, min(m,n))
+    elseif d == :U return triu(A.factors[1:min(m,n),1:n])
+    elseif d == :p
         p = [1:m]
-        for i in 1:length(A.ipiv)
+        for i in 1:length(A.pivots)
             tmp = p[i]
-            p[i] = p[A.ipiv[i]]
-            p[A.ipiv[i]] = tmp
+            p[i] = p[A.pivots[i]]
+            p[A.pivots[i]] = tmp
         end
         return p
-    end
-    if d == :P
+    elseif d == :P
         p = A[:p]
         P = zeros(T, m, m)
         for i in 1:m
@@ -186,13 +168,13 @@ function det{T}(A::LU{T})
     m, n = size(A)
     if m != n throw(DimensionMismatch("Matrix must be square")) end
     if A.info > 0; return zero(typeof(A.factors[1])); end
-    prod(diag(A.factors)) * (bool(sum(A.ipiv .!= 1:n) % 2) ? -one(T) : one(T))
+    prod(diag(A.factors)) * (bool(sum(A.pivots .!= 1:n) % 2) ? -one(T) : one(T))
 end
 
 function logdet2{T<:Real}(A::LU{T})  # return log(abs(det)) and sign(det)
     m, n = size(A); if m!=n error("matrix must be square") end
     dg = diag(A.factors)
-    s = (bool(sum(A.ipiv .!= 1:n) % 2) ? -one(T) : one(T)) * prod(sign(dg))
+    s = (bool(sum(A.pivots .!= 1:n) % 2) ? -one(T) : one(T)) * prod(sign(dg))
     return sum(log(abs(dg))) , s 
 end
 
@@ -204,7 +186,7 @@ end
 
 function logdet{T<:Complex}(A::LU{T})
     m, n = size(A); if m!=n error("matrix must be square") end
-    s = sum(log(diag(A.factors))) + (bool(sum(A.ipiv .!= 1:n) % 2) ? complex(0,pi) : 0) 
+    s = sum(log(diag(A.factors))) + (bool(sum(A.pivots .!= 1:n) % 2) ? complex(0,pi) : 0) 
     r,a = reim(s); a = a % 2pi; if a>pi a -=2pi elseif a<=-pi a+=2pi end
     return complex(r,a)    
 end
@@ -212,34 +194,34 @@ end
 
 function (\){T<:BlasFloat}(A::LU{T}, B::StridedVecOrMat{T})
     if A.info > 0; throw(SingularException(A.info)); end
-    LAPACK.getrs!('N', A.factors, A.ipiv, copy(B))
+    LAPACK.getrs!('N', A.factors, A.pivots, copy(B))
 end
 
 function At_ldiv_B{T<:BlasFloat}(A::LU{T}, B::StridedVecOrMat{T})
     if A.info > 0; throw(SingularException(A.info)); end
-    LAPACK.getrs!('T', A.factors, A.ipiv, copy(B))
+    LAPACK.getrs!('T', A.factors, A.pivots, copy(B))
 end
 
 function Ac_ldiv_B{T<:BlasComplex}(A::LU{T}, B::StridedVecOrMat{T})
     if A.info > 0; throw(SingularException(A.info)); end
-    LAPACK.getrs!('C', A.factors, A.ipiv, copy(B))
+    LAPACK.getrs!('C', A.factors, A.pivots, copy(B))
 end
 
 function At_ldiv_Bt{T<:BlasFloat}(A::LU{T}, B::StridedVecOrMat{T})
     if A.info > 0; throw(SingularException(A.info)); end
-    LAPACK.getrs!('T', A.factors, A.ipiv, transpose(B))
+    LAPACK.getrs!('T', A.factors, A.pivots, transpose(B))
 end
 
 function Ac_ldiv_Bc{T<:BlasComplex}(A::LU{T}, B::StridedVecOrMat{T})
     if A.info > 0; throw(SingularException(A.info)); end
-    LAPACK.getrs!('C', A.factors, A.ipiv, ctranspose(B))
+    LAPACK.getrs!('C', A.factors, A.pivots, ctranspose(B))
 end
 
 (/){T}(B::Matrix{T},A::LU{T}) = At_ldiv_Bt(A,B).'
 
 function inv(A::LU)
     if A.info > 0; return throw(SingularException(A.info)); end
-    LAPACK.getri!(copy(A.factors), A.ipiv)
+    LAPACK.getri!(copy(A.factors), A.pivots)
 end
 
 cond(A::LU, p) = 1.0/LinAlg.LAPACK.gecon!(p == 1 ? '1' : 'I', A.factors, norm(A[:L][A[:p],:]*A[:U], p))
@@ -298,14 +280,7 @@ function getindex{T<:BlasFloat}(A::QRPivoted{T}, d::Symbol)
     if d == :R return triu(A.hh[1:minimum(size(A)),:])
     elseif d == :Q return QRPivotedQ(A)
     elseif d == :p return A.pivots
-    elseif d == :P #The permutation matrix for the pivot
-        p = A[:p]
-        n = length(p)
-        P = zeros(T, n, n)
-        for i in 1:n
-            P[p[i],i] = one(T)
-        end
-        return P
+    elseif d == :P return perm_matrix(A.pivots, T) #The permutation matrix for the pivot
     end
     error("No such type field")
 end

From 18d0ec9a9e1493d19ef670d8a305e4b6e1e8bbf6 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jiahao=20Chen=20=28=E9=99=88=E5=AE=B6=E8=B1=AA=29?=
 <jiahao@mit.edu>
Date: Fri, 25 Oct 2013 21:36:47 -0400
Subject: [PATCH 3/8] Fix dimension finding to work for StridedMatOrVec

Revert A.hh[1] < eps to A.hh[1] == 0 testing for (\)(QRPivoted)
---
 base/linalg/factorization.jl | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/base/linalg/factorization.jl b/base/linalg/factorization.jl
index 2a95dda5548c0..0352d4283ef51 100644
--- a/base/linalg/factorization.jl
+++ b/base/linalg/factorization.jl
@@ -295,7 +295,7 @@ print_matrix(io::IO, A::Union(QRPackedQ, QRPivotedQ), rows::Integer, cols::Integ
 
 ## Multiplication by Q from the QR decomposition
 function *{T<:BlasFloat}(A::QRPackedQ{T}, B::StridedVecOrMat{T})
-    m, n = size(B)
+    m, n = size(B, 1), ndims(B)==1 ? 1 : size(B, 2)
     if m == size(A.vs, 1)
         Bc = copy(B)
     elseif m == size(A.vs, 2)
@@ -307,7 +307,7 @@ function *{T<:BlasFloat}(A::QRPackedQ{T}, B::StridedVecOrMat{T})
 end
 *{T<:BlasFloat}(A::StridedVecOrMat{T}, B::QRPackedQ{T}) = LAPACK.gemqrt!('R', 'N', B.vs, B.T, copy(A))
 function *{T<:BlasFloat}(A::QRPivotedQ{T}, B::StridedVecOrMat{T})
-    m, n = size(B)
+    m, n = size(B, 1), ndims(B)==1 ? 1 : size(B, 2)
     if m == size(A.hh, 1)
         Bc = copy(B)
     elseif m == size(A.hh, 2)
@@ -324,7 +324,7 @@ Ac_mul_B{T<:BlasComplex}(A::QRPackedQ{T},  B::StridedVecOrMat) = LAPACK.gemqrt!(
 Ac_mul_B{T<:BlasReal   }(A::QRPivotedQ{T}, B::StridedVecOrMat) = LAPACK.ormqr! ('L','T',A.hh,A.tau,copy(B))
 Ac_mul_B{T<:BlasComplex}(A::QRPivotedQ{T}, B::StridedVecOrMat) = LAPACK.ormqr! ('L','C',A.hh,A.tau,copy(B))
 function A_mul_Bc{T<:BlasFloat}(A::StridedVecOrMat{T}, B::QRPackedQ{T})
-    m, n = size(A)
+    m, n = size(A, 1), ndims(A)==1 ? 1 : size(A, 2)
     if n == size(B.vs, 1)
         Ac = copy(A)
     elseif n == size(B.vs, 2)
@@ -335,7 +335,7 @@ function A_mul_Bc{T<:BlasFloat}(A::StridedVecOrMat{T}, B::QRPackedQ{T})
     LAPACK.gemqrt!('R', iseltype(B.vs,Complex) ? 'C' : 'T', B.vs, B.T, Ac)
 end
 function A_mul_Bc{T<:BlasFloat}(A::StridedVecOrMat{T}, B::QRPivotedQ{T})
-    m, n = size(A)
+    m, n = size(A, 1), ndims(A)==1 ? 1 : size(A, 2)
     if n == size(B.hh, 1)
         Ac = copy(A)
     elseif n == size(B.hh, 2)
@@ -355,7 +355,7 @@ function (\){T<:BlasFloat}(A::QRPivoted{T}, B::StridedMatrix{T}, rcond::Real)
     nr, nrhs = minimum(size(A.hh)), size(B, 2)
     if nr == 0 return zeros(0, nrhs), 0 end
     ar = abs(A.hh[1])
-    if ar < eps(typeof(ar)) return zeros(nr, nrhs), 0 end
+    if ar == 0 return zeros(nr, nrhs), 0 end #XXX is 0 the correct underflow threshold?
     rnk = 1
     xmin, xmax = ones(T, nr), ones(T, nr)
     tmin, tmax = ar, ar

From 028542a93e9803ac1c291bfe8f31c6e6825cfb71 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jiahao=20Chen=20=28=E9=99=88=E5=AE=B6=E8=B1=AA=29?=
 <jiahao@mit.edu>
Date: Wed, 23 Oct 2013 19:20:22 -0400
Subject: [PATCH 4/8] Refactors QR code

- Deprecates qrp, qrpfact and qrpfact!
- Adds trailing pivoted::Bool option to qr, qrfact and qrpfact!
- Reorganizes QR factorizations section so that methods common to QR,
  QRPackedQ, QRPivoted, QRPivotedQ are grouped together.
- Remove improperly cleaned up blob from merge
- Replace methods defining default parameters with methods with default
  parameter values
- Rename QRPivoted.jpvt to QRPivoted.pivots
---
 base/deprecated.jl           |   4 +
 base/linalg/dense.jl         |   4 +-
 base/linalg/factorization.jl | 206 +++++++++++++++++++----------------
 3 files changed, 116 insertions(+), 98 deletions(-)

diff --git a/base/deprecated.jl b/base/deprecated.jl
index 89311493d1ed1..925131bfb0529 100644
--- a/base/deprecated.jl
+++ b/base/deprecated.jl
@@ -168,6 +168,10 @@ export PipeString
 #   copy!(dest::AbstractArray, doffs::Integer, src::Integer)
 # in abstractarray.jl
 @deprecate copy!(dest::AbstractArray, src, doffs::Integer)  copy!(dest, doffs, src)
+@deprecate qrpfact!(A)         qrfact!(A, true)
+@deprecate qrpfact(A)          qrfact(A, true)
+@deprecate qrp(A, thin)        qr(A, thin, true)
+@deprecate qrp(A)              qr(A, true)
 
 deprecated_ls() = run(`ls -l`)
 deprecated_ls(args::Cmd) = run(`ls -l $args`)
diff --git a/base/linalg/dense.jl b/base/linalg/dense.jl
index e1e78d55a3037..4ccbc25a32648 100644
--- a/base/linalg/dense.jl
+++ b/base/linalg/dense.jl
@@ -421,7 +421,7 @@ function factorize!{T}(A::Matrix{T})
         end
         return lufact!(A)
     end
-    return qrpfact!(A)
+    return qrfact!(A, true)
 end
 
 factorize(A::AbstractMatrix) = factorize!(copy(A))
@@ -441,7 +441,7 @@ function (\){T<:BlasFloat}(A::StridedMatrix{T}, B::StridedVecOrMat{T})
         istriu(A) && return \(Triangular(A, :U),B)
         return \(lufact(A),B)
     end
-    return qrpfact(A)\B
+    return qrfact(A, true)\B
 end
 
 ## Moore-Penrose inverse
diff --git a/base/linalg/factorization.jl b/base/linalg/factorization.jl
index ae1903be6409d..a17fee45cdef7 100644
--- a/base/linalg/factorization.jl
+++ b/base/linalg/factorization.jl
@@ -217,85 +217,60 @@ inv(A::LU)=@assertnonsingular LAPACK.getri!(copy(A.factors), A.ipiv) A.info
 
 cond(A::LU, p) = 1.0/LinAlg.LAPACK.gecon!(p == 1 ? '1' : 'I', A.factors, norm(A[:L][A[:p],:]*A[:U], p))
 
-## QR decomposition without column pivots. By the faster geqrt3
+####################
+# QR factorization #
+####################
+
 type QR{S<:BlasFloat} <: Factorization{S}
     vs::Matrix{S}                     # the elements on and above the diagonal contain the N-by-N upper triangular matrix R; the elements below the diagonal are the columns of V
     T::Matrix{S}                      # upper triangular factor of the block reflector.
 end
-QR{T<:BlasFloat}(A::StridedMatrix{T}, nb::Integer = min(minimum(size(A)), 36)) = QR(LAPACK.geqrt!(A, nb)...)
-
-qrfact!{T<:BlasFloat}(A::StridedMatrix{T}, args::Integer...) = QR(A, args...)
-qrfact!(A::StridedMatrix, args::Integer...) = qrfact!(float(A), args...)
-qrfact{T<:BlasFloat}(A::StridedMatrix{T}, args::Integer...) = qrfact!(copy(A), args...)
-qrfact(A::StridedMatrix, args::Integer...) = qrfact!(float(A), args...)
-qrfact(x::Integer) = qrfact(float(x))
-qrfact(x::Number) = QR(fill(one(x), 1, 1), fill(x, 1, 1))
-
-function qr(A::Union(Number, AbstractMatrix), thin::Bool=true)
-    F = qrfact(A)
-    full(F[:Q], thin), F[:R]
-end
 
-size(A::QR, args::Integer...) = size(A.vs, args...)
-
-function getindex(A::QR, d::Symbol)
-    d == :R && return triu(A.vs[1:minimum(size(A)),:])
-    d == :Q && return QRPackedQ(A)
-    throw(KeyError(d))
-end
+# QR decomposition without column pivots. By the faster geqrt3
+QR{T<:BlasFloat}(A::StridedMatrix{T}, nb::Integer = min(minimum(size(A)), 36)) = QR(LAPACK.geqrt!(A, nb)...)
 
 type QRPackedQ{S} <: AbstractMatrix{S} 
-    vs::Matrix{S}                      
-    T::Matrix{S}                       
+    vs::Matrix{S}
+    T::Matrix{S}
 end
 QRPackedQ(A::QR) = QRPackedQ(A.vs, A.T)
 
-size(A::QRPackedQ, dim::Integer) = 0 < dim ? (dim <= 2 ? size(A.vs, 1) : 1) : throw(BoundsError())
-size(A::QRPackedQ) = size(A, 1), size(A, 2)
-
-full{T<:BlasFloat}(A::QRPackedQ{T}, thin::Bool=true) = A * (thin ? eye(T, size(A.vs)...) : eye(T, size(A.vs,1)))
-
-print_matrix(io::IO, A::QRPackedQ, rows::Integer, cols::Integer) = print_matrix(io, full(A, false), rows, cols)
-
-## Multiplication by Q from the QR decomposition
-function *{T<:BlasFloat}(A::QRPackedQ{T}, B::StridedVecOrMat{T})
-    Bpad = size(B, 1)==size(A.vs, 2) ? [B; zeros(T, size(A.vs, 1) - size(A.vs, 2), size(B, 2))] : copy(B)
-    LAPACK.gemqrt!('L', 'N', A.vs, A.T, Bpad)
-end
-*{T<:BlasFloat}(A::StridedVecOrMat{T}, B::QRPackedQ{T}) = LAPACK.gemqrt!('R', 'N', B.vs, B.T, copy(A))
-Ac_mul_B{T<:BlasReal}(A::QRPackedQ{T}, B::StridedVecOrMat{T}) = LAPACK.gemqrt!('L','T',A.vs,A.T,copy(B))
-Ac_mul_B{T<:BlasComplex}(A::QRPackedQ{T}, B::StridedVecOrMat{T}) = LAPACK.gemqrt!('L','C',A.vs,A.T,copy(B))
-Ac_mul_B(A::QRPackedQ, B::StridedVecOrMat) = Ac_mul_B(A, float(B))
-function A_mul_Bc{T<:BlasFloat}(A::StridedVecOrMat{T}, B::QRPackedQ{T})
-    Apad = size(A, 2)==size(B.vs, 2) ? [A zeros(T, size(A, 1), size(B.vs, 1) - size(B.vs, 2))] : copy(A)     
-    LAPACK.gemqrt!('R', iseltype(B.vs,Complex) ? 'C' : 'T', B.vs, B.T, Apad)
-end
-## Least squares solution.  Should be more careful about cases with m < n
-\(A::QR, B::StridedVector) = Triangular(A[:R], :U)\(A[:Q]'B)[1:size(A, 2)]
-\(A::QR, B::StridedMatrix) = Triangular(A[:R], :U)\(A[:Q]'B)[1:size(A, 2),:]
-
 type QRPivoted{T} <: Factorization{T}
     hh::Matrix{T}
     tau::Vector{T}
-    jpvt::Vector{BlasInt}
+    pivots::Vector{BlasInt}
 end
 
-qrpfact!{T<:BlasFloat}(A::StridedMatrix{T}) = QRPivoted{T}(LAPACK.geqp3!(A)...)
-qrpfact!(A::StridedMatrix) = qrpfact!(float(A))
-qrpfact{T<:BlasFloat}(A::StridedMatrix{T}) = qrpfact!(copy(A))
-qrpfact(A::StridedMatrix) = qrpfact!(float(A))
+type QRPivotedQ{T} <: AbstractMatrix{T}
+    hh::Matrix{T}                      # Householder transformations and R
+    tau::Vector{T}                     # Scalar factors of transformations
+end
+QRPivotedQ(A::QRPivoted) = QRPivotedQ(A.hh, A.tau)
 
-function qrp(A::AbstractMatrix, thin::Bool=false)
-    F = qrpfact(A)
-    full(F[:Q], thin), F[:R], F[:p]
+function qr(A::Union(Number, AbstractMatrix), thin::Bool=true, pivoted::Bool=false)
+  F = qrfact(A, pivoted)
+  return pivoted ? (full(F[:Q], thin), F[:R], F[:p]) : (full(F[:Q], thin), F[:R])
 end
 
-size(A::QRPivoted, args::Integer...) = size(A.hh, args...)
+qrfact!{T<:BlasFloat}(A::StridedMatrix{T}, pivoted::Bool=false, args::Integer...) = pivoted ? QRPivoted{T}(LAPACK.geqp3!(A)...) : QR(A, args...) 
+qrfact!(A::StridedMatrix, pivoted::Bool=false) = qrfact!(float(A), pivoted)
+qrfact{T<:BlasFloat}(A::StridedMatrix{T}, pivoted::Bool=false, args::Integer...) = qrfact!(copy(A), pivoted, args...)
+qrfact(A::StridedMatrix, pivoted::Bool=false, args::Integer...) = qrfact!(float(A), pivoted, args...)
+qrfact(x::Integer, pivoted::Bool=false) = qrfact(float(x), pivoted)
+qrfact(x::Number, pivoted::Bool=false) = pivoted ? error("Not implemented") : QR(fill(one(x), 1, 1), fill(x, 1, 1))
+size(A::Union(QR,        QRPackedQ ), args::Integer...) = size(A.vs, args...)
+size(A::Union(QRPivoted, QRPivotedQ), args::Integer...) = size(A.hh, args...)
+
+function getindex(A::QR, d::Symbol)
+    d == :R && return triu(A.vs[1:minimum(size(A)),:])
+    d == :Q && return QRPackedQ(A)
+    throw(KeyError(d))
+end
 
 function getindex{T<:BlasFloat}(A::QRPivoted{T}, d::Symbol)
     d == :R && return triu(A.hh[1:minimum(size(A)),:])
     d == :Q && return QRPivotedQ(A)
-    d == :p && return A.jpvt
+    d == :p && return A.pivots
     if d == :P
         p = A[:p]
         n = length(p)
@@ -308,28 +283,94 @@ function getindex{T<:BlasFloat}(A::QRPivoted{T}, d::Symbol)
     throw(KeyError(d))
 end
 
+full{T<:BlasFloat}(A::QRPackedQ{T}, thin::Bool) = A * eye(T, thin ? size(A.T, 2) : size(A, 1))
+function full{T<:BlasFloat}(A::QRPivotedQ{T}, thin::Bool=true)
+    m, n = size(A.hh)
+    return LAPACK.orgqr!(thin ? copy(A.hh) : [A.hh zeros(T, m, max(0, m-n))], A.tau)
+end
+
+print_matrix(io::IO, A::Union(QRPackedQ, QRPivotedQ), rows::Integer, cols::Integer) = print_matrix(io, full(A), rows, cols)
+
+## Multiplication by Q from the QR decomposition
+function *{T<:BlasFloat}(A::QRPackedQ{T}, B::StridedVecOrMat{T})
+    m, n = size(B)
+    if m == size(A.vs, 1)
+        Bc = copy(B)
+    elseif m == size(A.vs, 2)
+        Bc = [B; zeros(T, size(A.vs,1)-m, n)]
+    else
+        throw(DimensionMismatch(""))
+    end
+    LAPACK.gemqrt!('L', 'N', A.vs, A.T, Bc)
+end
+*{T<:BlasFloat}(A::StridedVecOrMat{T}, B::QRPackedQ{T}) = LAPACK.gemqrt!('R', 'N', B.vs, B.T, copy(A))
+function *{T<:BlasFloat}(A::QRPivotedQ{T}, B::StridedVecOrMat{T})
+    m, n = size(B)
+    if m == size(A.hh, 1)
+        Bc = copy(B)
+    elseif m == size(A.hh, 2)
+        Bc = [B; zeros(T, size(A.hh,1)-m, n)]
+    else
+        throw(DimensionMismatch(""))
+    end
+    LAPACK.ormqr!('L', 'N', A.hh, A.tau, Bc)
+end
+*(A::StridedVecOrMat, B::QRPivotedQ) = LAPACK.ormqr!('R', 'N', B.hh, B.tau, copy(A))
+
+Ac_mul_B{T<:BlasReal   }(A::QRPackedQ{T},  B::StridedVecOrMat) = LAPACK.gemqrt!('L','T',A.vs,A.T,  copy(B))
+Ac_mul_B{T<:BlasComplex}(A::QRPackedQ{T},  B::StridedVecOrMat) = LAPACK.gemqrt!('L','C',A.vs,A.T,  copy(B))
+Ac_mul_B{T<:BlasReal   }(A::QRPivotedQ{T}, B::StridedVecOrMat) = LAPACK.ormqr! ('L','T',A.hh,A.tau,copy(B))
+Ac_mul_B{T<:BlasComplex}(A::QRPivotedQ{T}, B::StridedVecOrMat) = LAPACK.ormqr! ('L','C',A.hh,A.tau,copy(B))
+function A_mul_Bc{T<:BlasFloat}(A::StridedVecOrMat{T}, B::QRPackedQ{T})
+    m, n = size(A)
+    if n == size(B.vs, 1)
+        Ac = copy(A)
+    elseif n == size(B.vs, 2)
+        Ac = [B zeros(T, m, size(B.vs, 1) - n)]
+    else
+        throw(DimensionMismatch(""))
+    end
+    LAPACK.gemqrt!('R', iseltype(B.vs,Complex) ? 'C' : 'T', B.vs, B.T, Ac)
+end
+function A_mul_Bc{T<:BlasFloat}(A::StridedVecOrMat{T}, B::QRPivotedQ{T})
+    m, n = size(A)
+    if n == size(B.hh, 1)
+        Ac = copy(A)
+    elseif n == size(B.hh, 2)
+        Ac = [B zeros(T, m, size(B.hh, 1) - n)]
+    else
+        throw(DimensionMismatch(""))
+    end
+    LAPACK.ormqr!('R', iseltype(B.hh,Complex) ? 'C' : 'T', B.hh, B.tau, Ac)
+end
+
+## Least squares solution.
+#XXX Should be more careful about cases with m < n
+(\)(A::QR, B::StridedVector) = Triangular(A[:R], :U)\(A[:Q]'B)[1:size(A, 2)]
+(\)(A::QR, B::StridedMatrix) = Triangular(A[:R], :U)\(A[:Q]'B)[1:size(A, 2),:]
 # Julia implementation similarly to xgelsy
-function \{T<:BlasFloat}(A::QRPivoted{T}, B::StridedMatrix{T}, rcond::Real)
+function (\){T<:BlasFloat}(A::QRPivoted{T}, B::StridedMatrix{T}, rcond::Real)
     nr = minimum(size(A.hh))
     nrhs = size(B, 2)
     if nr == 0 return zeros(0, nrhs), 0 end
     ar = abs(A.hh[1])
-    if ar == 0 return zeros(nr, nrhs), 0 end
+    if ar < eps(typeof(ar)) return zeros(nr, nrhs), 0 end
     rnk = 1
     xmin = ones(T, nr)
     xmax = ones(T, nr)
     tmin = tmax = ar
     while rnk < nr
-        tmin, smin, cmin = LAPACK.laic1!(2, sub(xmin, 1:rnk), tmin, sub(A.hh, 1:rnk, rnk + 1), A.hh[rnk + 1, rnk + 1])
-        tmax, smax, cmax = LAPACK.laic1!(1, sub(xmax, 1:rnk), tmax, sub(A.hh, 1:rnk, rnk + 1), A.hh[rnk + 1, rnk + 1])
-        if tmax*rcond > tmin && break end
-        xmin[1:rnk + 1] = [smin*sub(xmin, 1:rnk), cmin]
-        xmax[1:rnk + 1] = [smax*sub(xmin, 1:rnk), cmax]
+        Ahh = sub(A.hh, 1:rnk, rnk+1), A.hh[rnk+1, rnk+1]
+        tmin, smin, cmin = LAPACK.laic1!(2, sub(xmin, 1:rnk), tmin, Ahh...)
+        tmax, smax, cmax = LAPACK.laic1!(1, sub(xmax, 1:rnk), tmax, Ahh...)
+        if tmax*rcond > tmin break end
+        xmin[1:rnk+1] = [smin*sub(xmin, 1:rnk), cmin]
+        xmax[1:rnk+1] = [smax*sub(xmin, 1:rnk), cmax] #XXX should this be xmax?
         rnk += 1
         # if cond(r[1:rnk, 1:rnk])*rcond < 1 break end
     end
     C, tau = LAPACK.tzrzf!(A.hh[1:rnk,:])
-    X = [Triangular(C[1:rnk,1:rnk],:U)\(A[:Q]'B)[1:rnk,:]; zeros(T, size(A.hh, 2) - rnk, nrhs)]
+    X = [Triangular(C[1:rnk,1:rnk],:U)\(A[:Q]'B)[1:rnk,:]; zeros(T, size(A.hh,2)-rnk, nrhs)]
     LAPACK.ormrz!('L', iseltype(B, Complex) ? 'C' : 'T', C, tau, X)
     return X[invperm(A[:p]),:], rnk
 end
@@ -337,41 +378,14 @@ end
 \(A::QRPivoted, B::StridedMatrix) = (\)(A, B, sqrt(eps(typeof(float(real(B[1]))))))[1]
 \(A::QRPivoted, B::StridedVector) = (\)(A, reshape(B, length(B), 1))[:]
 
-type QRPivotedQ{T} <: AbstractMatrix{T}
-    hh::Matrix{T}                       # Householder transformations and R
-    tau::Vector{T}                      # Scalar factors of transformations
-end
-QRPivotedQ(A::QRPivoted) = QRPivotedQ(A.hh, A.tau)
-
-size(A::QRPivotedQ, dims::Integer) = dims > 0 ? (dims < 3 ? size(A.hh, 1) : 1) : throw(BoundsError())
-
-function full{T<:BlasFloat}(A::QRPivotedQ{T}, thin::Bool=true)
-    m, n = size(A.hh)
-    Ahhpad = thin ? copy(A.hh) : [A.hh zeros(T, m, max(0, m - n))]
-    LAPACK.orgqr!(Ahhpad, A.tau)
-end
-
-print_matrix(io::IO, A::QRPivotedQ, rows::Integer, cols::Integer) = print_matrix(io, full(A, false), rows, cols)
-
-## Multiplication by Q from the Pivoted QR decomposition
-function *{T<:BlasFloat}(A::QRPivotedQ{T}, B::StridedVecOrMat{T})
-    Bpad = size(A.hh, 2)==size(B, 1) ? [B; zeros(T, size(A.hh, 1) - size(A.hh, 2), size(B, 2))] : copy(B)
-    LAPACK.ormqr!('L', 'N', A.hh, A.tau, Bpad)
-end
-
-Ac_mul_B{T<:BlasReal}(A::QRPivotedQ{T}, B::StridedVecOrMat{T}) = LAPACK.ormqr!('L','T',A.hh,A.tau,copy(B))
-Ac_mul_B{T<:BlasComplex}(A::QRPivotedQ{T}, B::StridedVecOrMat{T}) = LAPACK.ormqr!('L','C',A.hh,A.tau,copy(B))
-Ac_mul_B(A::QRPivotedQ, B::StridedVecOrMat) = Ac_mul_B(A, float(B))
-*(A::StridedVecOrMat, B::QRPivotedQ) = LAPACK.ormqr!('R', 'N', B.hh, B.tau, copy(A))
-function A_mul_Bc{T<:BlasFloat}(A::StridedVecOrMat{T}, B::QRPivotedQ{T})
-    Apad = size(A, 2)==size(B.hh, 2) ? [A zeros(T, size(A, 1), size(B.hh, 1) - size(B.hh, 2))] : copy(A)
-    LAPACK.ormqr!('R', iseltype(B.hh,Complex) ? 'C' : 'T', B.hh, B.tau, Apad)
-end
-
 ##TODO:  Add methods for rank(A::QRP{T}) and adjust the (\) method accordingly
 ##       Add rcond methods for Cholesky, LU, QR and QRP types
 ## Lower priority: Add LQ, QL and RQ factorizations
 
+############################
+# Hessenberg factorization #
+############################
+
 # FIXME! Should add balancing option through xgebal
 type Hessenberg{T} <: Factorization{T}
     hh::Matrix{T}

From f30f4a7b456f93049e7669077828e7616db13f75 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jiahao=20Chen=20=28=E9=99=88=E5=AE=B6=E8=B1=AA=29?=
 <jiahao@mit.edu>
Date: Thu, 24 Oct 2013 17:38:36 -0400
Subject: [PATCH 5/8] Deprecates cholpfact, cholpfact!

---
 base/deprecated.jl           |   4 +
 base/linalg/factorization.jl | 185 +++++++++++++++++------------------
 2 files changed, 95 insertions(+), 94 deletions(-)

diff --git a/base/deprecated.jl b/base/deprecated.jl
index 925131bfb0529..6d80f529465e3 100644
--- a/base/deprecated.jl
+++ b/base/deprecated.jl
@@ -172,6 +172,10 @@ export PipeString
 @deprecate qrpfact(A)          qrfact(A, true)
 @deprecate qrp(A, thin)        qr(A, thin, true)
 @deprecate qrp(A)              qr(A, true)
+@deprecate cholpfact!(A)            cholfact!(A, :U, true)
+@deprecate cholpfact!(A,tol)        cholfact!(A, :U, true, tol)
+@deprecate cholpfact!(A,uplo,tol)   cholfact!(A, uplo, true, tol)
+@deprecate cholpfact(A)             cholfact(A, :U, true)
 
 deprecated_ls() = run(`ls -l`)
 deprecated_ls(args::Cmd) = run(`ls -l $args`)
diff --git a/base/linalg/factorization.jl b/base/linalg/factorization.jl
index a17fee45cdef7..612d095d8b4c0 100644
--- a/base/linalg/factorization.jl
+++ b/base/linalg/factorization.jl
@@ -10,6 +10,16 @@ macro assertnonsingular(A, info)
    :(($info)==0 ? $A : throw(SingularException($info)))
 end
 
+# Generate permutation matrix from a permutation
+function perm_matrix{T<:Integer}(A::Vector{T})
+    n = length(A)
+    P = zeros(T, n, n)
+    for i=1:n
+        P[A[i], i] = one(output_type)
+    end
+    P
+end
+
 \(F::Factorization, B::AbstractVecOrMat) = A_ldiv_B!(F, copy(B))
 Ac_ldiv_B(F::Factorization, B::AbstractVecOrMat) = Ac_ldiv_B!(F, copy(B))
 At_ldiv_B(F::Factorization, B::AbstractVecOrMat) = At_ldiv_B!(F, copy(B))
@@ -17,27 +27,47 @@ A_ldiv_B!(F::Factorization, B::StridedVecOrMat) = A_ldiv_B!(F, float(B))
 Ac_ldiv_B!(F::Factorization, B::StridedVecOrMat) = Ac_ldiv_B!(F, float(B))
 At_ldiv_B!(F::Factorization, B::StridedVecOrMat) = At_ldiv_B!(F, float(B))
 
+##########################
+# Cholesky factorization #
+##########################
+
 type Cholesky{T<:BlasFloat} <: Factorization{T}
     UL::Matrix{T}
     uplo::Char
 end
 
-function cholfact!{T<:BlasFloat}(A::StridedMatrix{T}, uplo::Symbol)
+type CholeskyPivoted{T<:BlasFloat} <: Factorization{T}
+    UL::Matrix{T}
+    uplo::Char
+    piv::Vector{BlasInt}
+    rank::BlasInt
+    tol::Real
+    info::BlasInt
+end
+
+function CholeskyPivoted{T<:BlasFloat}(A::StridedMatrix{T}, uplo::Char, tol::Real)
+    A, piv, rank, info = LAPACK.pstrf!(uplo, A, tol)
+    CholeskyPivoted{T}(uplo=='U' ? triu!(A) : tril!(A), uplo, piv, rank, tol, info)
+end
+
+function cholfact!{T<:BlasFloat}(A::StridedMatrix{T}, uplo::Symbol=:U, pivoted::Bool=false,
+      tol::Real=-1.0)
+  if pivoted
+    CholeskyPivoted(A, string(uplo)[1], tol)
+  else
     uplochar = string(uplo)[1]
     C, info = LAPACK.potrf!(uplochar, A)
-    @assertposdef Cholesky(C, uplochar) info
+    info==0 ? Cholesky(C, uplochar) : throw(PosDefException(info))
+  end
 end
-cholfact!(A::StridedMatrix, args...) = cholfact!(float(A), args...)
-cholfact!{T<:BlasFloat}(A::StridedMatrix{T}) = cholfact!(A, :U)
-cholfact{T<:BlasFloat}(A::StridedMatrix{T}, args...) = cholfact!(copy(A), args...)
-cholfact(A::StridedMatrix, args...) = cholfact!(float(A), args...)
-cholfact(x::Number) = @assertposdef Cholesky(fill(sqrt(x), 1, 1), :U) !(imag(x) == 0 && real(x) > 0)
 
-chol(A::Union(Number, AbstractMatrix), uplo::Symbol) = cholfact(A, uplo)[uplo]
-chol(A::Union(Number, AbstractMatrix)) = cholfact(A, :U)[:U]
+cholfact!(A::StridedMatrix, uplo::Symbol=:U, pivoted::Bool=false, args...) = cholfact!(float(A), uplo, pivoted, args...)
+cholfact{T<:BlasFloat}(A::StridedMatrix{T}, uplo::Symbol=:U, pivoted::Bool=false, args...) = cholfact!(copy(A), uplo, pivoted, args...)
+cholfact(A::StridedMatrix, uplo::Symbol=:U, pivoted::Bool=false, args...) = cholfact!(float(A), uplo, pivoted, args...)
+cholfact(A::StridedMatrix, pivoted::Bool=false, args...) = cholfact(float(A), :U, pivoted, args...)
+cholfact(x::Number, pivoted::Bool=false) = imag(x) == 0 && real(x) > 0 ? Cholesky(fill(sqrt(x), 1, 1), :U) : throw(PosDefException(1))
 
-size(C::Cholesky) = size(C.UL)
-size(C::Cholesky,d::Integer) = size(C.UL,d)
+chol(A::Union(Number, AbstractMatrix), uplo::Symbol=:U) = cholfact(A, uplo)[uplo]
 
 function getindex(C::Cholesky, d::Symbol)
     C.uplo == 'U' ? triu!(C.UL) : tril!(C.UL)
@@ -49,70 +79,28 @@ function getindex(C::Cholesky, d::Symbol)
     throw(KeyError(d))
 end
 
-A_ldiv_B!{T<:BlasFloat}(C::Cholesky{T}, B::StridedVecOrMat{T}) = LAPACK.potrs!(C.uplo, C.UL, B)
-
-function det{T}(C::Cholesky{T})
-    dd = one(T)
-    for i in 1:size(C.UL,1) dd *= abs2(C.UL[i,i]) end
-    dd
-end
-
-function logdet{T}(C::Cholesky{T})
-    dd = zero(T)
-    for i in 1:size(C.UL,1) dd += log(C.UL[i,i]) end
-    dd + dd # instead of 2.0dd which can change the type
-end
-
-inv(C::Cholesky)=symmetrize_conj!(LAPACK.potri!(C.uplo, copy(C.UL)), C.uplo)
-
-## Pivoted Cholesky
-type CholeskyPivoted{T<:BlasFloat} <: Factorization{T}
-    UL::Matrix{T}
-    uplo::Char
-    piv::Vector{BlasInt}
-    rank::BlasInt
-    tol::Real
-    info::BlasInt
-end
-function CholeskyPivoted{T<:BlasFloat}(A::StridedMatrix{T}, uplo::Char, tol::Real)
-    A, piv, rank, info = LAPACK.pstrf!(uplo, A, tol)
-    CholeskyPivoted{T}((uplo == 'U' ? triu! : tril!)(A), uplo, piv, rank, tol, info)
-end
-
-chkfullrank(C::CholeskyPivoted) = C.rank<size(C.UL, 1) && throw(RankDeficientException(C.info))
-
-cholpfact!(A::StridedMatrix, args...) = cholpfact!(float(A), args...)
-cholpfact!{T<:BlasFloat}(A::StridedMatrix{T}, uplo::Symbol, tol::Real) = CholeskyPivoted(A, string(uplo)[1], tol)
-cholpfact!{T<:BlasFloat}(A::StridedMatrix{T}, tol::Real) = cholpfact!(A, :U, tol)
-cholpfact!{T<:BlasFloat}(A::StridedMatrix{T}) = cholpfact!(A, -1.)
-cholpfact{T<:BlasFloat}(A::StridedMatrix{T}, args...) = cholpfact!(copy(A), args...)
-cholpfact(A::StridedMatrix, args...) = cholpfact!(float(A), args...)
-
-size(C::CholeskyPivoted) = size(C.UL)
-size(C::CholeskyPivoted,d::Integer) = size(C.UL,d)
-
 getindex(C::CholeskyPivoted) = C.UL, C.piv
+
 function getindex{T<:BlasFloat}(C::CholeskyPivoted{T}, d::Symbol)
     (d == :U || d == :L) && return symbol(C.uplo) == d ? C.UL : C.UL'
     d == :p && return C.piv
-    if d == :P
-        n = size(C, 1)
-        P = zeros(T, n, n)
-        for i=1:n
-            P[C.piv[i],i] = one(T)
-        end
-        return P
-    end
+    d == :P && return perm_matrix(C.piv)
     throw(KeyError(d))
 end
 
+size(C::Union(Cholesky,CholeskyPivoted)) = size(C.UL)
+size(C::Union(Cholesky,CholeskyPivoted),d::Integer) = size(C.UL,d)
+
+rank(C::CholeskyPivoted) = C.rank
+rank(C::Cholesky) = size(C.UL, 1) #full rank
+
+A_ldiv_B!{T<:BlasFloat}(C::Cholesky{T}, B::StridedVecOrMat{T}) = LAPACK.potrs!(C.uplo, C.UL, B)
 function A_ldiv_B!{T<:BlasFloat}(C::CholeskyPivoted{T}, B::StridedVector{T})
-    chkfullrank(C)
+    C.rank < size(C.UL, 1) ? throw(RankDeficientException(C.info)) :
     ipermute!(LAPACK.potrs!(C.uplo, C.UL, permute!(B, C.piv)), C.piv)
 end
-
 function A_ldiv_B!{T<:BlasFloat}(C::CholeskyPivoted{T}, B::StridedMatrix{T})
-    chkfullrank(C)
+    C.rank < size(C.UL, 1) ? throw(RankDeficientException(C.info)) : nothing
     n = size(C, 1)
     for i=1:size(B, 2)
         permute!(sub(B, 1:n, i), C.piv)
@@ -124,9 +112,16 @@ function A_ldiv_B!{T<:BlasFloat}(C::CholeskyPivoted{T}, B::StridedMatrix{T})
     B
 end
 
-rank(C::CholeskyPivoted) = C.rank
+det{T}(C::Cholesky{T}) = prod(abs2(diag(C.UL)))
+det{T}(C::CholeskyPivoted{T}) = C.rank < size(C.UL, 1) ? real(zero(T)) : prod(abs2(diag(C.UL)))
+
+logdet{T}(C::Cholesky{T}) = 2sum(log(diag(C.UL)))
+logdet{T}(C::CholeskyPivoted{T}) = C.rank < size(C.UL, 1) ? real(convert(T, -Inf)) : 2sum(log(diag(C.UL)))
 
-det{T}(C::CholeskyPivoted{T}) = C.rank<size(C.UL,1) ? real(zero(T)) : prod(abs2(diag(C.UL)))
+function inv(C::Cholesky)
+    Ci, info = LAPACK.potri!(C.uplo, copy(C.UL))
+    info == 0 ? symmetrize_conj!(Ci, C.uplo) : throw(SingularException(info))
+end
 
 function inv(C::CholeskyPivoted)
     chkfullrank(C)
@@ -134,13 +129,20 @@ function inv(C::CholeskyPivoted)
     (symmetrize!(LAPACK.potri!(C.uplo, copy(C.UL)), C.uplo))[ipiv, ipiv]
 end
 
-## LU
+####################
+# LU factorization #
+####################
+
 type LU{T<:BlasFloat} <: Factorization{T}
     factors::Matrix{T}
-    ipiv::Vector{BlasInt}
+    pivots::Vector{BlasInt}
     info::BlasInt
 end
-LU{T<:BlasFloat}(A::StridedMatrix{T}) = LU{T}(LAPACK.getrf!(A)...)
+
+function LU{T<:BlasFloat}(A::StridedMatrix{T})
+    factors, pivots, info = LAPACK.getrf!(A)
+    LU{T}(factors, pivots, info)
+end
 
 lufact!(A::StridedMatrix) = lufact!(float(A))
 lufact!{T<:BlasFloat}(A::StridedMatrix{T}) = LU(A)
@@ -178,16 +180,17 @@ function getindex{T}(A::LU{T}, d::Symbol)
 end
 
 function det{T}(A::LU{T})
-    n = chksquare(A)
+    m, n = size(A)
+    m != n && throw(DimensionMismatch("Matrix must be square"))
     A.info > 0 && return zero(typeof(A.factors[1]))
-    return prod(diag(A.factors)) * (bool(sum(A.ipiv .!= 1:n) % 2) ? -one(T) : one(T))
+    prod(diag(A.factors)) * (bool(sum(A.pivots .!= 1:n) % 2) ? -one(T) : one(T))
 end
 
 function logdet2{T<:Real}(A::LU{T})  # return log(abs(det)) and sign(det)
     n = chksquare(A)
     dg = diag(A.factors)
-    s = (bool(sum(A.ipiv .!= 1:n) % 2) ? -one(T) : one(T)) * prod(sign(dg))
-    sum(log(abs(dg))), s 
+    s = (bool(sum(A.pivots .!= 1:n) % 2) ? -one(T) : one(T)) * prod(sign(dg))
+    return sum(log(abs(dg))), s 
 end
 
 function logdet{T<:Real}(A::LU{T})
@@ -197,23 +200,25 @@ function logdet{T<:Real}(A::LU{T})
 end
 
 function logdet{T<:Complex}(A::LU{T})
-    n = chksquare(A)
-    s = sum(log(diag(A.factors))) + (bool(sum(A.ipiv .!= 1:n) % 2) ? complex(0,pi) : 0) 
-    r, a = reim(s)
-    a = pi-mod(pi-a,2pi) #Take principal branch with argument (-pi,pi] 
-    complex(r,a)    
+    m, n = size(A)
+    m != n && throw(DimensionMisMatch("matrix must be square"))
+    s = sum(log(diag(A.factors))) + (bool(sum(A.pivots .!= 1:n) % 2) ? complex(0,pi) : 0) 
+    r, a = reim(s); a = a % 2pi; if a>pi a -=2pi elseif a<=-pi a+=2pi end
+    return complex(r,a)    
 end
 
-
-A_ldiv_B!{T<:BlasFloat}(A::LU{T}, B::StridedVecOrMat{T}) = @assertnonsingular LAPACK.getrs!('N', A.factors, A.ipiv, B) A.info
-At_ldiv_B{T<:BlasFloat}(A::LU{T}, B::StridedVecOrMat{T}) = @assertnonsingular LAPACK.getrs!('T', A.factors, A.ipiv, copy(B)) A.info
-Ac_ldiv_B{T<:BlasComplex}(A::LU{T}, B::StridedVecOrMat{T}) = @assertnonsingular LAPACK.getrs!('C', A.factors, A.ipiv, copy(B)) A.info
-At_ldiv_Bt{T<:BlasFloat}(A::LU{T}, B::StridedVecOrMat{T}) = @assertnonsingular LAPACK.getrs!('T', A.factors, A.ipiv, transpose(B)) A.info
-Ac_ldiv_Bc{T<:BlasComplex}(A::LU{T}, B::StridedVecOrMat{T}) = @assertnonsingular LAPACK.getrs!('C', A.factors, A.ipiv, ctranspose(B)) A.info
+A_ldiv_B!{T<:BlasFloat}(A::LU{T}, B::StridedVecOrMat{T}) = @assertnonsingular LAPACK.getrs!('N', A.factors, A.pivots, B) A.info
+At_ldiv_B{T<:BlasFloat}(A::LU{T}, B::StridedVecOrMat{T}) = @assertnonsingular LAPACK.getrs!('T', A.factors, A.pivots, copy(B)) A.info
+Ac_ldiv_B{T<:BlasComplex}(A::LU{T}, B::StridedVecOrMat{T}) = @assertnonsingular LAPACK.getrs!('C', A.factors, A.pivots, copy(B)) A.info
+At_ldiv_Bt{T<:BlasFloat}(A::LU{T}, B::StridedVecOrMat{T}) = @assertnonsingular LAPACK.getrs!('T', A.factors, A.pivots, transpose(B)) A.info
+Ac_ldiv_Bc{T<:BlasComplex}(A::LU{T}, B::StridedVecOrMat{T}) = @assertnonsingular LAPACK.getrs!('C', A.factors, A.pivots, ctranspose(B)) A.info
 
 /{T}(B::Matrix{T},A::LU{T}) = At_ldiv_Bt(A,B).'
 
-inv(A::LU)=@assertnonsingular LAPACK.getri!(copy(A.factors), A.ipiv) A.info
+function inv(A::LU)
+    if A.info > 0; return throw(SingularException(A.info)); end
+    LAPACK.getri!(copy(A.factors), A.pivots)
+end
 
 cond(A::LU, p) = 1.0/LinAlg.LAPACK.gecon!(p == 1 ? '1' : 'I', A.factors, norm(A[:L][A[:p],:]*A[:U], p))
 
@@ -271,15 +276,7 @@ function getindex{T<:BlasFloat}(A::QRPivoted{T}, d::Symbol)
     d == :R && return triu(A.hh[1:minimum(size(A)),:])
     d == :Q && return QRPivotedQ(A)
     d == :p && return A.pivots
-    if d == :P
-        p = A[:p]
-        n = length(p)
-        P = zeros(T, n, n)
-        for i in 1:n
-            P[p[i],i] = one(T)
-        end
-        return P
-    end
+    d == :P && return perm_matrix(A.pivots) #The permutation matrix for the pivot
     throw(KeyError(d))
 end
 

From 34d54f59a351e291e60df24ac238dcc224c7975c Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jiahao=20Chen=20=28=E9=99=88=E5=AE=B6=E8=B1=AA=29?=
 <jiahao@mit.edu>
Date: Fri, 25 Oct 2013 21:36:47 -0400
Subject: [PATCH 6/8] Fix dimension finding to work for StridedMatOrVec

Revert A.hh[1] < eps to A.hh[1] == 0 testing for (\)(QRPivoted)
---
 base/linalg/factorization.jl | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/base/linalg/factorization.jl b/base/linalg/factorization.jl
index 612d095d8b4c0..3f6b76f40395b 100644
--- a/base/linalg/factorization.jl
+++ b/base/linalg/factorization.jl
@@ -290,7 +290,7 @@ print_matrix(io::IO, A::Union(QRPackedQ, QRPivotedQ), rows::Integer, cols::Integ
 
 ## Multiplication by Q from the QR decomposition
 function *{T<:BlasFloat}(A::QRPackedQ{T}, B::StridedVecOrMat{T})
-    m, n = size(B)
+    m, n = size(B, 1), ndims(B)==1 ? 1 : size(B, 2)
     if m == size(A.vs, 1)
         Bc = copy(B)
     elseif m == size(A.vs, 2)
@@ -302,7 +302,7 @@ function *{T<:BlasFloat}(A::QRPackedQ{T}, B::StridedVecOrMat{T})
 end
 *{T<:BlasFloat}(A::StridedVecOrMat{T}, B::QRPackedQ{T}) = LAPACK.gemqrt!('R', 'N', B.vs, B.T, copy(A))
 function *{T<:BlasFloat}(A::QRPivotedQ{T}, B::StridedVecOrMat{T})
-    m, n = size(B)
+    m, n = size(B, 1), ndims(B)==1 ? 1 : size(B, 2)
     if m == size(A.hh, 1)
         Bc = copy(B)
     elseif m == size(A.hh, 2)
@@ -319,7 +319,7 @@ Ac_mul_B{T<:BlasComplex}(A::QRPackedQ{T},  B::StridedVecOrMat) = LAPACK.gemqrt!(
 Ac_mul_B{T<:BlasReal   }(A::QRPivotedQ{T}, B::StridedVecOrMat) = LAPACK.ormqr! ('L','T',A.hh,A.tau,copy(B))
 Ac_mul_B{T<:BlasComplex}(A::QRPivotedQ{T}, B::StridedVecOrMat) = LAPACK.ormqr! ('L','C',A.hh,A.tau,copy(B))
 function A_mul_Bc{T<:BlasFloat}(A::StridedVecOrMat{T}, B::QRPackedQ{T})
-    m, n = size(A)
+    m, n = size(A, 1), ndims(A)==1 ? 1 : size(A, 2)
     if n == size(B.vs, 1)
         Ac = copy(A)
     elseif n == size(B.vs, 2)
@@ -330,7 +330,7 @@ function A_mul_Bc{T<:BlasFloat}(A::StridedVecOrMat{T}, B::QRPackedQ{T})
     LAPACK.gemqrt!('R', iseltype(B.vs,Complex) ? 'C' : 'T', B.vs, B.T, Ac)
 end
 function A_mul_Bc{T<:BlasFloat}(A::StridedVecOrMat{T}, B::QRPivotedQ{T})
-    m, n = size(A)
+    m, n = size(A, 1), ndims(A)==1 ? 1 : size(A, 2)
     if n == size(B.hh, 1)
         Ac = copy(A)
     elseif n == size(B.hh, 2)
@@ -351,7 +351,7 @@ function (\){T<:BlasFloat}(A::QRPivoted{T}, B::StridedMatrix{T}, rcond::Real)
     nrhs = size(B, 2)
     if nr == 0 return zeros(0, nrhs), 0 end
     ar = abs(A.hh[1])
-    if ar < eps(typeof(ar)) return zeros(nr, nrhs), 0 end
+    if ar == 0 return zeros(nr, nrhs), 0 end #XXX is 0 the correct underflow threshold?
     rnk = 1
     xmin = ones(T, nr)
     xmax = ones(T, nr)

From 2a9e42da78b43bbb67b6b514da9e8cd87fa53b8e Mon Sep 17 00:00:00 2001
From: "Viral B. Shah" <viral@mayin.org>
Date: Sun, 29 Dec 2013 11:22:22 +0530
Subject: [PATCH 7/8] WIP

---
 base/linalg.jl               | 1 +
 base/linalg/factorization.jl | 4 ++--
 base/linalg/matmul.jl        | 4 ++++
 3 files changed, 7 insertions(+), 2 deletions(-)

diff --git a/base/linalg.jl b/base/linalg.jl
index 8968705571b7d..89cd6da2c60ed 100644
--- a/base/linalg.jl
+++ b/base/linalg.jl
@@ -111,6 +111,7 @@ export
     svdvals!,
     svdvals,
     symmetrize!,
+    symmetrize_conj!,
     trace,
     transpose,
     tril,
diff --git a/base/linalg/factorization.jl b/base/linalg/factorization.jl
index 3f6b76f40395b..47bb2d63a1afd 100644
--- a/base/linalg/factorization.jl
+++ b/base/linalg/factorization.jl
@@ -164,8 +164,8 @@ function getindex{T}(A::LU{T}, d::Symbol)
     d == :U && return triu(A.factors[1:min(m,n),1:n])
     if d == :p
         p = [1:m]
-        for i in 1:length(A.ipiv)
-            p[i], p[A.ipiv[i]] = p[A.ipiv[i]], p[i]
+        for i in 1:length(A.pivots)
+            p[i], p[A.pivots[i]] = p[A.pivots[i]], p[i]
         end
         return p
     elseif d == :P
diff --git a/base/linalg/matmul.jl b/base/linalg/matmul.jl
index 47f6d4aaeedfb..78f728c762b47 100644
--- a/base/linalg/matmul.jl
+++ b/base/linalg/matmul.jl
@@ -147,6 +147,8 @@ function symmetrize!(A::StridedMatrix, uplo::Char='U')
     A
 end
 
+symmetrize!(A::Number, uplo) = A
+
 function symmetrize_conj!(A::StridedMatrix, uplo::Char='U')
     n = chksquare(A)
     @chkuplo
@@ -162,6 +164,8 @@ function symmetrize_conj!(A::StridedMatrix, uplo::Char='U')
     A
 end
 
+symmetrize_conj!(A::Number, uplo) = A
+
 function gemv{T<:BlasFloat}(y::StridedVector{T}, tA::Char, A::StridedMatrix{T}, x::StridedVector{T})
     stride(A, 1)==1 || return generic_matvecmul(y, tA, A, x)
 

From 5fed276c07872106c681350028516314fd581fc1 Mon Sep 17 00:00:00 2001
From: Andreas Noack Jensen <andreasnoackjensen@gmail.com>
Date: Mon, 30 Dec 2013 11:26:43 +0100
Subject: [PATCH 8/8] Add chkfullrank. Fix inv(Cholesky). Define promotion
 rules for Q multiplication. Update tests.

---
 base/linalg/factorization.jl | 27 +++++++++++++++++++++------
 test/linalg.jl               |  2 +-
 2 files changed, 22 insertions(+), 7 deletions(-)

diff --git a/base/linalg/factorization.jl b/base/linalg/factorization.jl
index cb93698d6ab6c..66590e588159f 100644
--- a/base/linalg/factorization.jl
+++ b/base/linalg/factorization.jl
@@ -119,10 +119,12 @@ logdet{T}(C::Cholesky{T}) = 2sum(log(diag(C.UL)))
 logdet{T}(C::CholeskyPivoted{T}) = C.rank < size(C.UL, 1) ? real(convert(T, -Inf)) : 2sum(log(diag(C.UL)))
 
 function inv(C::Cholesky)
-    Ci, info = LAPACK.potri!(C.uplo, copy(C.UL))
-    info == 0 ? symmetrize_conj!(Ci, C.uplo) : throw(SingularException(info))
+    Ci = LAPACK.potri!(C.uplo, copy(C.UL))
+    symmetrize_conj!(Ci, C.uplo)
 end
 
+chkfullrank(C::CholeskyPivoted) = C.rank == size(C,1) ? true : throw(SingularException(C.rank))
+
 function inv(C::CholeskyPivoted)
     chkfullrank(C)
     ipiv = invperm(C.piv)
@@ -314,10 +316,10 @@ function *{T<:BlasFloat}(A::QRPivotedQ{T}, B::StridedVecOrMat{T})
 end
 *(A::StridedVecOrMat, B::QRPivotedQ) = LAPACK.ormqr!('R', 'N', B.hh, B.tau, copy(A))
 
-Ac_mul_B{T<:BlasReal   }(A::QRPackedQ{T},  B::StridedVecOrMat) = LAPACK.gemqrt!('L','T',A.vs,A.T,  copy(B))
-Ac_mul_B{T<:BlasComplex}(A::QRPackedQ{T},  B::StridedVecOrMat) = LAPACK.gemqrt!('L','C',A.vs,A.T,  copy(B))
-Ac_mul_B{T<:BlasReal   }(A::QRPivotedQ{T}, B::StridedVecOrMat) = LAPACK.ormqr! ('L','T',A.hh,A.tau,copy(B))
-Ac_mul_B{T<:BlasComplex}(A::QRPivotedQ{T}, B::StridedVecOrMat) = LAPACK.ormqr! ('L','C',A.hh,A.tau,copy(B))
+Ac_mul_B{T<:BlasReal   }(A::QRPackedQ{T},  B::StridedVecOrMat{T}) = LAPACK.gemqrt!('L','T',A.vs,A.T,  copy(B))
+Ac_mul_B{T<:BlasComplex}(A::QRPackedQ{T},  B::StridedVecOrMat{T}) = LAPACK.gemqrt!('L','C',A.vs,A.T,  copy(B))
+Ac_mul_B{T<:BlasReal   }(A::QRPivotedQ{T}, B::StridedVecOrMat{T}) = LAPACK.ormqr! ('L','T',A.hh,A.tau,copy(B))
+Ac_mul_B{T<:BlasComplex}(A::QRPivotedQ{T}, B::StridedVecOrMat{T}) = LAPACK.ormqr! ('L','C',A.hh,A.tau,copy(B))
 function A_mul_Bc{T<:BlasFloat}(A::StridedVecOrMat{T}, B::QRPackedQ{T})
     m, n = size(A, 1), ndims(A)==1 ? 1 : size(A, 2)
     if n == size(B.vs, 1)
@@ -340,6 +342,19 @@ function A_mul_Bc{T<:BlasFloat}(A::StridedVecOrMat{T}, B::QRPivotedQ{T})
     end
     LAPACK.ormqr!('R', iseltype(B.hh,Complex) ? 'C' : 'T', B.hh, B.tau, Ac)
 end
+# promotions
+*{T<:BlasFloat}(A::QRPackedQ{T}, B::StridedVector) = *(A,convert(Vector{T},B))
+*{T<:BlasFloat}(A::QRPivotedQ{T}, B::StridedVector) = *(A,convert(Vector{T},B))
+*{T<:BlasFloat}(A::QRPackedQ{T}, B::StridedMatrix) = *(A,convert(Matrix{T},B))
+*{T<:BlasFloat}(A::QRPivotedQ{T}, B::StridedMatrix) = *(A,convert(Matrix{T},B))
+*{T<:BlasFloat}(A::StridedMatrix, B::QRPackedQ{T}) = *(convert(Matrix{T},A),B)
+*{T<:BlasFloat}(A::StridedMatrix, B::QRPivotedQ{T}) = *(convert(Matrix{T},A),B)
+Ac_mul_B{T<:BlasFloat}(A::QRPackedQ{T}, B::StridedVector) = Ac_mul_B(A,convert(Vector{T},B))
+Ac_mul_B{T<:BlasFloat}(A::QRPivotedQ{T}, B::StridedVector) = Ac_mul_B(A,convert(Vector{T},B))
+Ac_mul_B{T<:BlasFloat}(A::QRPackedQ{T}, B::StridedMatrix) = Ac_mul_B(A,convert(Matrix{T},B))
+Ac_mul_B{T<:BlasFloat}(A::QRPivotedQ{T}, B::StridedMatrix) = Ac_mul_B(A,convert(Matrix{T},B))
+A_mul_Bc{T<:BlasFloat}(A::StridedMatrix, B::QRPackedQ{T}) = Ac_mul_B(convert(Matrix{T},A),B)
+A_mul_Bc{T<:BlasFloat}(A::StridedMatrix, B::QRPivotedQ{T}) = Ac_mul_B(convert(Matrix{T},A),B)
 
 ## Least squares solution.
 #XXX Should be more careful about cases with m < n
diff --git a/test/linalg.jl b/test/linalg.jl
index 56b18659976cc..73847781ea969 100644
--- a/test/linalg.jl
+++ b/test/linalg.jl
@@ -49,7 +49,7 @@ for elty in (Float32, Float64, Complex64, Complex128, Int)
     @test_approx_eq l*l' apd
 
     # pivoted Choleksy decomposition
-    cpapd = cholpfact(apd)
+    cpapd = cholfact(apd, true)
     @test rank(cpapd) == n
     @test all(diff(diag(real(cpapd.UL))).<=0.) # diagonal should be non-increasing
     @test_approx_eq b apd * (cpapd\b)