Merge pull request #8403 from JuliaLang/teh/cholesky

Use constants for BLAS chars (speed up Cholesky factorization of small matrices)

Author: timholy

base/linalg/bunchkaufman.jl

@@ -11,12 +11,12 @@ end
 
 function bkfact!{T<:BlasReal}(A::StridedMatrix{T}, uplo::Symbol=:U, symmetric::Bool=issym(A))
     symmetric || error("The Bunch-Kaufman decomposition is only valid for symmetric matrices")
-    LD, ipiv = LAPACK.sytrf!(string(uplo)[1] , A)
-    BunchKaufman(LD, ipiv, string(uplo)[1], symmetric)
+    LD, ipiv = LAPACK.sytrf!(char_uplo(uplo) , A)
+    BunchKaufman(LD, ipiv, char_uplo(uplo), symmetric)
 end
 function bkfact!{T<:BlasComplex}(A::StridedMatrix{T}, uplo::Symbol=:U, symmetric::Bool=issym(A))
-    LD, ipiv = (symmetric ? LAPACK.sytrf! : LAPACK.hetrf!)(string(uplo)[1] , A)
-    BunchKaufman(LD, ipiv, string(uplo)[1], symmetric)
+    LD, ipiv = (symmetric ? LAPACK.sytrf! : LAPACK.hetrf!)(char_uplo(uplo) , A)
+    BunchKaufman(LD, ipiv, char_uplo(uplo), symmetric)
 end
 bkfact{T<:BlasFloat}(A::StridedMatrix{T}, uplo::Symbol=:U, symmetric::Bool=issym(A)) = bkfact!(copy(A), uplo, symmetric)
 bkfact{T}(A::StridedMatrix{T}, uplo::Symbol=:U, symmetric::Bool=issym(A)) = bkfact!(convert(Matrix{promote_type(Float32,typeof(sqrt(one(T))))},A),uplo,symmetric)

base/linalg/cholesky.jl

@@ -14,7 +14,7 @@ immutable CholeskyPivoted{T} <: Factorization{T}
 end
 
 function chol!{T<:BlasFloat}(A::StridedMatrix{T}, uplo::Symbol=:U)
-    C, info = LAPACK.potrf!(string(uplo)[1], A)
+    C, info = LAPACK.potrf!(char_uplo(uplo), A)
     return @assertposdef Triangular(C, uplo, false) info
 end
 
@@ -57,7 +57,7 @@ function chol!{T}(A::AbstractMatrix{T}, uplo::Symbol=:U)
 end
 
 function cholfact!{T<:BlasFloat}(A::StridedMatrix{T}, uplo::Symbol=:U; pivot=false, tol=0.0)
-    uplochar = string(uplo)[1]
+    uplochar = char_uplo(uplo)
     if pivot
         A, piv, rank, info = LAPACK.pstrf!(uplochar, A, tol)
         return CholeskyPivoted{T}(A, uplochar, piv, rank, tol, info)
@@ -66,6 +66,12 @@ function cholfact!{T<:BlasFloat}(A::StridedMatrix{T}, uplo::Symbol=:U; pivot=fal
 end
 cholfact!(A::AbstractMatrix, uplo::Symbol=:U) = Cholesky{eltype(A),typeof(A),uplo}(chol!(A, uplo).data)
 
+function cholfact!{T<:BlasFloat,S,UpLo}(C::Cholesky{T,S,UpLo})
+    _, info = LAPACK.potrf!(char_uplo(UpLo), C.UL)
+    info[1]>0 && throw(PosDefException(info[1]))
+    C
+end
+
 cholfact{T<:BlasFloat}(A::StridedMatrix{T}, uplo::Symbol=:U; pivot=false, tol=0.0) = cholfact!(copy(A), uplo, pivot=pivot, tol=tol)
 function cholfact{T}(A::StridedMatrix{T}, uplo::Symbol=:U; pivot=false, tol=0.0)
     S = promote_type(typeof(chol(one(T))),Float32)

base/linalg/dense.jl

@@ -20,7 +20,7 @@ scale!{T<:BlasFloat}(X::Array{T}, s::Number) = scale!(X, convert(T, s))
 scale!{T<:BlasComplex}(X::Array{T}, s::Real) = BLAS.scal!(length(X), oftype(real(zero(T)),s), X, 1)
 
 #Test whether a matrix is positive-definite
-isposdef!{T<:BlasFloat}(A::StridedMatrix{T}, UL::Symbol) = LAPACK.potrf!(string(UL)[1], A)[2] == 0
+isposdef!{T<:BlasFloat}(A::StridedMatrix{T}, UL::Symbol) = LAPACK.potrf!(char_uplo(UL), A)[2] == 0
 isposdef!(A::StridedMatrix) = ishermitian(A) && isposdef!(A, :U)
 
 isposdef{T}(A::AbstractMatrix{T}, UL::Symbol) = (S = typeof(sqrt(one(T))); isposdef!(S == T ? copy(A) : convert(AbstractMatrix{S}, A), UL))

base/linalg/symmetric.jl

@@ -3,12 +3,12 @@ immutable Symmetric{T,S<:AbstractMatrix{T}} <: AbstractMatrix{T}
     data::S
     uplo::Char
 end
-Symmetric(A::AbstractMatrix, uplo::Symbol=:U) = (chksquare(A);Symmetric{eltype(A),typeof(A)}(A, string(uplo)[1]))
+Symmetric(A::AbstractMatrix, uplo::Symbol=:U) = (chksquare(A);Symmetric{eltype(A),typeof(A)}(A, char_uplo(uplo)))
 immutable Hermitian{T,S<:AbstractMatrix{T}} <: AbstractMatrix{T}
     data::S
     uplo::Char
 end
-Hermitian(A::AbstractMatrix, uplo::Symbol=:U) = (chksquare(A);Hermitian{eltype(A),typeof(A)}(A, string(uplo)[1]))
+Hermitian(A::AbstractMatrix, uplo::Symbol=:U) = (chksquare(A);Hermitian{eltype(A),typeof(A)}(A, char_uplo(uplo)))
 typealias HermOrSym{T,S} Union(Hermitian{T,S}, Symmetric{T,S})
 typealias RealHermSymComplexHerm{T<:Real,S} Union(Hermitian{T,S}, Symmetric{T,S}, Hermitian{Complex{T},S})
 

base/linalg/triangular.jl

@@ -7,6 +7,10 @@ function Triangular{T}(A::AbstractMatrix{T}, uplo::Symbol, isunit::Bool=false)
     return Triangular{T,typeof(A),uplo,isunit}(A)
 end
 
+const CHARU = 'U'
+const CHARL = 'L'
+char_uplo(uplo::Symbol) = uplo == :U ? CHARU : (uplo == :L ? CHARL : throw(ArgumentError("uplo argument must be either :U or :L")))
+
 +{T, MT, uplo}(A::Triangular{T, MT, uplo, false}, B::Triangular{T, MT, uplo, false}) = Triangular(A.data + B.data, uplo)
 +{T, MT}(A::Triangular{T, MT, :U, false}, B::Triangular{T, MT, :U, true}) = Triangular(A.data + triu(B.data, 1) + I, :U)
 +{T, MT}(A::Triangular{T, MT, :L, false}, B::Triangular{T, MT, :L, true}) = Triangular(A.data + tril(B.data, -1) + I, :L)