Skip to content

Commit

Permalink
Reinstate lufact, cholfact, cholpfact, qrfact, and qrpfact.
Browse files Browse the repository at this point in the history 
Matlab-like factorizations are now provided in lu, chol, and qr.
  • Loading branch information
@ViralBShah
ViralBShah committed Mar 13, 2013
1 parent dbf37cc commit 96cd587
Show file tree
Hide file tree
Showing 6 changed files with 52 additions and 28 deletions.
6 changes: 0 additions & 6 deletions base/deprecated.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -146,12 +146,6 @@ end
@deprecate expr(hd, a...) Expr(hd, a...)
@deprecate expr(hd, a::Array{Any,1}) Expr(hd, a...)

@deprecate cholfact chol
@deprecate cholpfact cholp
@deprecate lufact lu
@deprecate qrfact qr
@deprecate qrpfact qrp

@deprecate logb exponent
@deprecate ref_shape index_shape
@deprecate assign_shape_check setindex_shape_check
Expand Down
4 changes: 4 additions & 0 deletions base/exports.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -550,7 +550,11 @@ export

# linear algebra
chol,
cholfact,
cholfact!,
cholp,
cholpfact,
cholpfact!,
cond,
cross,
ctranspose,
Expand Down
2 changes: 1 addition & 1 deletion base/linalg/dense.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -415,7 +415,7 @@ function det(A::Matrix)
end
det(x::Number) = x

logdet(A::Matrix) = 2.0 * sum(log(diag(chol(A)[:U])))
logdet(A::Matrix) = 2.0 * sum(log(diag(cholfact(A)[:U])))

function inv(A::StridedMatrix)
if istriu(A) return inv(Triangular(A, 'U')) end
Expand Down
52 changes: 37 additions & 15 deletions base/linalg/factorization.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -13,10 +13,14 @@ type CholeskyDense{T<:BlasFloat} <: Factorization{T}
end
CholeskyDense{T<:BlasFloat}(A::Matrix{T}, uplo::Char) = CholeskyDense{T}(A, uplo)

chol(A::Matrix, uplo::Symbol) = CholeskyDense(copy(A), string(uplo)[1])
chol(A::Matrix) = chol(A, :U)
chol{T<:Integer}(A::Matrix{T}, args...) = chol(float64(A), args...)
chol(x::Number) = imag(x) == 0 && real(x) > 0 ? sqrt(x) : error("Argument not positive-definite")
cholfact!(A::Matrix, uplo::Symbol) = CholeskyDense(A, string(uplo)[1])
cholfact(A::Matrix, uplo::Symbol) = cholfact!(copy(A), uplo)
cholfact!(A::Matrix) = cholfact!(A, :U)
cholfact(A::Matrix) = cholfact(A, :U)
cholfact{T<:Integer}(A::Matrix{T}, args...) = cholfact(float(A), args...)
cholfact(x::Number) = imag(x) == 0 && real(x) > 0 ? sqrt(x) : error("Argument not positive-definite")

chol(A) = cholfact(A, :U)[:U]

size(C::CholeskyDense) = size(C.UL)
size(C::CholeskyDense,d::Integer) = size(C.UL,d)
Expand Down Expand Up @@ -59,10 +63,13 @@ function CholeskyPivotedDense{T<:BlasFloat}(A::Matrix{T}, uplo::Char, tol::Real)
CholeskyPivotedDense{T}(uplo == 'U' ? triu!(A) : tril!(A), uplo, piv, rank, tol, info)
end

cholp(A::Matrix, uplo::Symbol, tol::Real) = CholeskyPivotedDense(copy(A), string(uplo)[1], tol)
cholp(A::Matrix, tol::Real) = cholp(A, :U, tol)
cholp(A::Matrix) = cholp(A, -1.)
cholp{T<:Int}(A::Matrix{T}, args...) = cholp(float64(A), args...)
cholpfact!(A::Matrix, uplo::Symbol, tol::Real) = CholeskyPivotedDense(A, string(uplo)[1], tol)
cholpfact(A::Matrix, uplo::Symbol, tol::Real) = cholpfact!(copy(A), uplo, tol)
cholpfact!(A::Matrix, tol::Real) = cholpfact!(A, :U, tol)
cholpfact(A::Matrix, tol::Real) = cholpfact(A, :U, tol)
cholpfact!(A::Matrix) = cholpfact!(A, -1.)
cholpfact(A::Matrix) = cholpfact(A, -1.)
cholpfact{T<:Int}(A::Matrix{T}, args...) = cholpfact(float(A), args...)

size(C::CholeskyPivotedDense) = size(C.UL)
size(C::CholeskyPivotedDense,d::Integer) = size(C.UL,d)
Expand Down Expand Up @@ -127,9 +134,16 @@ function LUDense{T<:BlasFloat}(A::Matrix{T})
LUDense{T}(LU, ipiv, info)
end

lu(A::Matrix) = LUDense(copy(A))
lu{T<:Integer}(A::Matrix{T}) = lu(float(A))
lu(x::Number) = (one(x), x, [1])
lufact!(A::Matrix) = LUDense(A)
lufact(A::Matrix) = lufact!(copy(A))
lufact!{T<:Integer}(A::Matrix{T}) = lufact!(float(A))
lufact{T<:Integer}(A::Matrix{T}) = lufact(float(A))
lufact(x::Number) = (one(x), x, [1])

function lu(A::Matrix)
F = lufact(A)
return (F[:L], F[:U], F[:P])
end

size(A::LUDense) = size(A.LU)
size(A::LUDense,n) = size(A.LU,n)
Expand Down Expand Up @@ -182,9 +196,15 @@ type QRDense{S} <: Factorization{S}
end
QRDense(A::Matrix) = QRDense(LAPACK.geqrt3!(A)...)

qr(A::Matrix) = QRDense(copy(A))
qr{T<:Integer}(A::Matrix{T}) = qr(float(A))
qr(x::Number) = (one(x), x)
qrfact!(A::Matrix) = QRDense(A)
qrfact(A::Matrix) = qrfact!(copy(A))
qrfact{T<:Integer}(A::Matrix{T}) = qrfact(float(A))
qrfact(x::Number) = (one(x), x)

function qr(A::Matrix)
F = qrfact(A)
return (F[:Q], F[:R])
end

size(A::QRDense, args::Integer...) = size(A.vs, args...)

Expand Down Expand Up @@ -254,7 +274,9 @@ type QRPivotedDense{T} <: Factorization{T}
end
end
QRPivotedDense{T<:BlasFloat}(A::Matrix{T}) = QRPivotedDense{T}(LAPACK.geqp3!(A)...)
qrp(A::Matrix) = QRPivotedDense(copy(A))

qrpfact!(A::Matrix) = QRPivotedDense(A)
qrpfact(A::Matrix) = qrpfact!(copy(A))
# QRDenseQ(A::QRPivotedDense) = QRDenseQ(A.hh, A.tau)

size(A::QRPivotedDense, args::Integer...) = size(A.hh, args...)
Expand Down
4 changes: 4 additions & 0 deletions base/linalg/linalg.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -26,7 +26,11 @@ export

# Functions
chol,
cholfact,
cholfact!,
cholp,
cholpfact,
cholpfact!,
cond,
copy!,
cross,
Expand Down
12 changes: 6 additions & 6 deletions test/linalg.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,18 +8,18 @@ for elty in (Float32, Float64, Complex64, Complex128)
apd = a'*a # symmetric positive-definite
b = convert(Vector{elty}, b)

capd = chol(apd) # upper Cholesky factor
capd = cholfact(apd) # upper Cholesky factor
r = capd[:U]
@test_approx_eq r'*r apd
@test_approx_eq b apd * (capd\b)
@test_approx_eq apd * inv(capd) eye(elty, n)
@test_approx_eq a*(capd\(a'*b)) b # least squares soln for square a
@test_approx_eq det(capd) det(apd)

l = chol(apd, :L)[:L] # lower Cholesky factor
l = cholfact(apd, :L)[:L] # lower Cholesky factor
@test_approx_eq l*l' apd

cpapd = cholp(apd) # pivoted Choleksy decomposition
cpapd = cholpfact(apd) # pivoted Choleksy decomposition
@test rank(cpapd) == n
@test all(diff(diag(real(cpapd.UL))).<=0.) # diagonal should be non-increasing
@test_approx_eq b apd * (cpapd\b)
Expand All @@ -32,21 +32,21 @@ for elty in (Float32, Float64, Complex64, Complex128)
@test_approx_eq inv(bc2) * apd eye(elty, n)
@test_approx_eq apd * (bc2\b) b

lua = lu(a) # LU decomposition
lua = lufact(a) # LU decomposition
l,u,p = lua[:L], lua[:U], lua[:p]
@test_approx_eq l*u a[p,:]
@test_approx_eq l[invperm(p),:]*u a
@test_approx_eq a * inv(lua) eye(elty, n)
@test_approx_eq a*(lua\b) b

qra = qr(a) # QR decomposition
qra = qrfact(a) # QR decomposition
q,r = qra[:Q], qra[:R]
@test_approx_eq q'*full(q, false) eye(elty, n)
@test_approx_eq q*full(q, false)' eye(elty, n)
@test_approx_eq q*r a
@test_approx_eq a*(qra\b) b

qrpa = qrp(a) # pivoted QR decomposition
qrpa = qrpfact(a) # pivoted QR decomposition
q,r,p = qrpa[:Q], qrpa[:R], qrpa[:p]
@test_approx_eq q'*full(q, false) eye(elty, n)
@test_approx_eq q*full(q, false)' eye(elty, n)
Expand Down

1 comment on commit 96cd587

@StefanKarpinski
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Man, what is going on here? Can you summarize how the decision was reached to do this and then revert it? I feel like this interface design is going in circles.

Please sign in to comment.