Merge pull request #5330 from JuliaLang/anj/pivots

Deprecate cholpfact and qrpfact and introduce keyword arguments for pivoting. Make thin a keyword argument.

NEWS.md

@@ -105,6 +105,9 @@ Deprecated or removed
     argument specifying the floating point type to which they apply. The old
     behaviour and `[get/set/with]_bigfloat_rounding` functions are deprecated ([#5007])
 
+  * cholpfact and qrpfact are deprecated in favor of keyword arguments in 
+    cholfact and qrfact.
+
 [#4042]: https://github.com/JuliaLang/julia/issues/4042
 [#5164]: https://github.com/JuliaLang/julia/issues/5164
 [#5214]: https://github.com/JuliaLang/julia/issues/5214
@@ -135,6 +138,7 @@ Deprecated or removed
 [#5277]: https://github.com/JuliaLang/julia/issues/5277
 [#987]: https://github.com/JuliaLang/julia/issues/987
 [#2345]: https://github.com/JuliaLang/julia/issues/2345
+[#5330]: https://github.com/JuliaLang/julia/issues/5330
 
 Julia v0.2.0 Release Notes
 ==========================

base/deprecated.jl

@@ -168,6 +168,14 @@ 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, pivot=true)
+@deprecate qrpfact(A)          qrfact(A, pivot=true)
+@deprecate qrp(A, thin)        qr(A, thin=thin, pivot=true)
+@deprecate qrp(A)              qr(A, pivot=true)
+@deprecate cholpfact!(A)            cholfact!(A, :U, pivot=true)
+@deprecate cholpfact!(A,tol=tol)        cholfact!(A, :U, pivot=true, tol=tol)
+@deprecate cholpfact!(A,uplo,tol=tol)   cholfact!(A, uplo, pivot=true, tol=tol)
+@deprecate cholpfact(A)             cholfact(A, :U, pivot=true)
 
 deprecated_ls() = run(`ls -l`)
 deprecated_ls(args::Cmd) = run(`ls -l $args`)

base/linalg.jl

@@ -111,6 +111,7 @@ export
     svdvals!,
     svdvals,
     symmetrize!,
+    symmetrize_conj!,
     trace,
     transpose,
     tril,

base/linalg/dense.jl

@@ -423,7 +423,7 @@ function factorize!{T}(A::Matrix{T})
         end
         return lufact!(A)
     end
-    return qrpfact!(A)
+    return qrfact!(A,pivot=true)
 end
 
 factorize(A::AbstractMatrix) = factorize!(copy(A))
@@ -443,7 +443,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,pivot=true)\B
 end
 
 ## Moore-Penrose inverse

base/linalg/factorization.jl

@@ -17,83 +17,51 @@ 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
+type CholeskyPivoted{T<:BlasFloat} <: Factorization{T}
+    UL::Matrix{T}
+    uplo::Char
+    piv::Vector{BlasInt}
+    rank::BlasInt
+    tol::Real
+    info::BlasInt
+end
 
-function cholfact!{T<:BlasFloat}(A::StridedMatrix{T}, uplo::Symbol)
+function cholfact!{T<:BlasFloat}(A::StridedMatrix{T}, uplo::Symbol=:U; pivot=false, tol=-1.0)
     uplochar = string(uplo)[1]
-    C, info = LAPACK.potrf!(uplochar, A)
-    @assertposdef Cholesky(C, uplochar) info
+    if pivot
+        A, piv, rank, info = LAPACK.pstrf!(uplochar, A, tol)
+        return CholeskyPivoted{T}(A, uplochar, piv, rank, tol, info)
+    else
+        C, info = LAPACK.potrf!(uplochar, A)
+        return @assertposdef Cholesky(C, uplochar) 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{T<:BlasFloat}(A::StridedMatrix{T}, uplo::Symbol=:U; pivot=false, tol=-1.0) = cholfact!(copy(A), uplo, pivot=pivot, tol=tol)
+cholfact(A::StridedMatrix, uplo::Symbol=:U; pivot=false, tol=-1.0) = cholfact!(float(A), uplo, pivot=pivot, tol=tol)
 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]
+chol(A::Union(Number, AbstractMatrix), uplo::Symbol=:U) = cholfact(A, uplo)[uplo]
 
 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
+    d == :U && return triu!(symbol(C.uplo) == d ? C.UL : C.UL')
+    d == :L && return tril!(symbol(C.uplo) == d ? C.UL : C.UL')
+    d == :UL && return Triangular(C.UL, C.uplo)
     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 == :U && return triu!(symbol(C.uplo) == d ? C.UL : C.UL')
+    d == :L && return tril!(symbol(C.uplo) == d ? C.UL : C.UL')
     d == :p && return C.piv
     if d == :P
         n = size(C, 1)
@@ -106,6 +74,10 @@ function getindex{T<:BlasFloat}(C::CholeskyPivoted{T}, d::Symbol)
     throw(KeyError(d))
 end
 
+show(io::IO, C::Cholesky) = (println("$(typeof(C)) with factor:");show(io,C[symbol(C.uplo)]))
+
+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)
     ipermute!(LAPACK.potrs!(C.uplo, C.UL, permute!(B, C.piv)), C.piv)
@@ -124,17 +96,35 @@ function A_ldiv_B!{T<:BlasFloat}(C::CholeskyPivoted{T}, B::StridedMatrix{T})
     B
 end
 
-rank(C::CholeskyPivoted) = C.rank
+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
 
 det{T}(C::CholeskyPivoted{T}) = C.rank<size(C.UL,1) ? real(zero(T)) : prod(abs2(diag(C.UL)))
 
+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)
+
 function inv(C::CholeskyPivoted)
     chkfullrank(C)
     ipiv = invperm(C.piv)
     (symmetrize!(LAPACK.potri!(C.uplo, copy(C.UL)), C.uplo))[ipiv, ipiv]
 end
 
-## LU
+chkfullrank(C::CholeskyPivoted) = C.rank<size(C.UL, 1) && throw(RankDeficientException(C.info))
+
+rank(C::CholeskyPivoted) = C.rank
+
+####################
+# LU Factorization #
+####################
 type LU{T<:BlasFloat} <: Factorization{T}
     factors::Matrix{T}
     ipiv::Vector{BlasInt}
@@ -204,7 +194,6 @@ function logdet{T<:Complex}(A::LU{T})
     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
@@ -217,23 +206,33 @@ inv(A::LU)=@assertnonsingular LAPACK.getri!(copy(A.factors), A.ipiv) A.info
 
 cond(A::LU, p) = inv(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 #
+####################
+
+# Note. For QR factorization without pivoting, the WY representation based method introduced in LAPACK 3.4
 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.
+    vs::Matrix{S}
+    T::Matrix{S}
 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...)
+type QRPivoted{T} <: Factorization{T}
+    hh::Matrix{T}
+    tau::Vector{T}
+    jpvt::Vector{BlasInt}
+end
+
+qrfact!{T<:BlasFloat}(A::StridedMatrix{T}; pivot=false) = pivot ? QRPivoted{T}(LAPACK.geqp3!(A)...) : QR(A)
+qrfact!(A::StridedMatrix; pivot=false) = qrfact!(float(A), pivot=pivot)
+qrfact{T<:BlasFloat}(A::StridedMatrix{T}; pivot=false) = qrfact!(copy(A), pivot=pivot)
+qrfact(A::StridedMatrix; pivot=false) = qrfact!(float(A), pivot=pivot)
 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]
+function qr(A::Union(Number, AbstractMatrix); pivot=false, thin::Bool=true)
+    F = qrfact(A, pivot=pivot)
+    full(F[:Q], thin=thin), F[:R]
 end
 
 size(A::QR, args::Integer...) = size(A.vs, args...)
@@ -253,9 +252,9 @@ 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)))
+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)
+print_matrix(io::IO, A::QRPackedQ, rows::Integer, cols::Integer) = print_matrix(io, full(A, thin=false), rows, cols)
 
 ## Multiplication by Q from the QR decomposition
 function *{T<:BlasFloat}(A::QRPackedQ{T}, B::StridedVecOrMat{T})
@@ -274,22 +273,6 @@ end
 \(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=false)
-    F = qrpfact(A)
-    full(F[:Q], thin), F[:R], F[:p]
-end
-
 size(A::QRPivoted, args::Integer...) = size(A.hh, args...)
 
 function getindex{T<:BlasFloat}(A::QRPivoted{T}, d::Symbol)
@@ -345,13 +328,13 @@ 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)
+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)
+print_matrix(io::IO, A::QRPivotedQ, rows::Integer, cols::Integer) = print_matrix(io, full(A, thin=false), rows, cols)
 
 ## Multiplication by Q from the Pivoted QR decomposition
 function *{T<:BlasFloat}(A::QRPivotedQ{T}, B::StridedVecOrMat{T})

base/linalg/lapack.jl

@@ -1674,7 +1674,7 @@ for (posv, potrf, potri, potrs, pstrf, elty, rtyp) in
                   (Ptr{BlasChar}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, Ptr{BlasInt},
                    Ptr{BlasInt}, Ptr{$rtyp}, Ptr{$rtyp}, Ptr{BlasInt}),
                   &uplo, &n, A, &max(1,stride(A,2)), piv, rank, &tol, work, info)
-            @lapackerror
+            @assertargsok
             A, piv, rank[1], info[1] #Stored in PivotedCholesky
         end
     end

base/linalg/triangular.jl

@@ -82,11 +82,9 @@ end
 
 size(A::Triangular, args...) = size(A.UL, args...)
 convert(::Type{Matrix}, A::Triangular) = full(A)
-full(A::Triangular) = A.UL
+full(A::Triangular) = A.uplo == 'U' ? triu!(A.UL) : tril!(A.UL)
 
-getindex{T}(A::Triangular{T}, i::Integer, j::Integer) = i == j ? A.UL[i,j] : ((A.uplo == 'U') == (i < j) ? getindex(A.UL, i, j) : zero(T))
-
-print_matrix(io::IO, A::Triangular, rows::Integer, cols::Integer) = print_matrix(io, full(A), rows, cols)
+getindex{T}(A::Triangular{T}, i::Integer, j::Integer) = i == j ? (A.unitdiag == 'U' ? one(T) : A.UL[i,j]) : ((A.uplo == 'U') == (i < j) ? getindex(A.UL, i, j) : zero(T))
 
 istril(A::Triangular) = A.uplo == 'L' || istriu(A.UL)
 istriu(A::Triangular) = A.uplo == 'U' || istril(A.UL)