diff --git a/base/export.jl b/base/export.jl index c81a82ca00f65..923a23e48b5af 100644 --- a/base/export.jl +++ b/base/export.jl @@ -90,13 +90,12 @@ export Zip, Stat, Factorization, - Cholesky, - LU, + CholeskyDense, + LUDense, LUTridiagonal, - LDLT, LDLTTridiagonal, - QR, - QRP, + QRDense, + QRPDense, # Exceptions ArgumentError, @@ -623,12 +622,14 @@ export eig, expm, eye, + factors, ishermitian, issym, issym_rnd, istril, istriu, kron, + ldlt, linreg, lu, lu!, @@ -636,9 +637,11 @@ export matmul3x3, norm, qr, + qrp, randsym, rank, rref, + sdd, solve, svd, svdvals, @@ -647,8 +650,8 @@ export trideig, tril, triu, - qrp, - sdd, + tril!, + triu!, # dequeues append!, diff --git a/base/factorizations.jl b/base/factorizations.jl index 6311c5e35fbd2..d655aad5d4724 100644 --- a/base/factorizations.jl +++ b/base/factorizations.jl @@ -112,50 +112,95 @@ end abstract Factorization{T} -type Cholesky{T} <: Factorization{T} +type CholeskyDense{T} <: Factorization{T} LR::Matrix{T} uplo::LapackChar - function Cholesky(A::Matrix{T}, ul::LapackChar) - if ul != 'U' && ul != 'L' error("Cholesky: uplo must be 'U' or 'L'") end - Acopy = copy(A) - _jl_lapack_potrf(ul, Acopy) == 0 ? new(ul == 'U' ? triu(Acopy) : tril(Acopy), ul) : error("Cholesky: Matrix is not positive-definite") +end + +# chol() does not check that input matrix is symmetric/hermitian +# It simply uses upper triangular half +function chol{T<:LapackScalar}(A::Matrix{T}, ul::LapackChar) + if ul != 'U' && ul != 'L'; error("Cholesky: uplo must be 'U' or 'L'"); end + C = CholeskyDense{T}(Array(T, size(A)), ul) + chol(C, A) +end +chol{T<:Real}(A::Matrix{T}, ul::LapackChar) = chol(float64(A), ul) +chol(A) = chol(A, 'U') + +function _chol{T<:LapackScalar}(C::CholeskyDense{T}) + if _jl_lapack_potrf(C.uplo, C.LR) != 0 + error("Cholesky: Matrix is not positive-definite") end + if C.uplo == 'U' + triu!(C.LR) + else + tril!(C.LR) + end + C +end +function chol{T}(C::CholeskyDense{T}, A::Matrix{T}) + copy_to(C.LR, A) + _chol(C) +end +function chol!{T<:LapackScalar}(A::Matrix{T}, ul::LapackChar) + C = CholeskyDense{T}(A, ul) + _chol(C) end + +factors(C::CholeskyDense) = C.LR -Cholesky(A::Matrix) = Cholesky{eltype(A)}(A, 'U') -(\){T<:Union(Float64,Float32,Complex128,Complex64)}(C::Cholesky{T}, B::StridedVecOrMat{T}) = +(\){T<:LapackScalar}(C::CholeskyDense{T}, B::StridedVecOrMat{T}) = _jl_lapack_potrs(C.uplo, C.LR, copy(B)) -inv(C::Cholesky) = _jl_lapack_potri(C.uplo, copy(C.LR)) # should symmetrize the result +inv{T<:LapackScalar}(C::CholeskyDense{T}) = _jl_lapack_potri(C.uplo, copy(C.LR)) # should symmetrize the result -type LU{T} <: Factorization{T} +type LUDense{T} <: Factorization{T} lu::Matrix{T} ipiv::Vector{Int32} - function LU(lu::Matrix{T}, ipiv::Vector{Int32}) + function LUDense(lu::Matrix{T}, ipiv::Vector{Int32}) m, n = size(lu) - m == numel(ipiv) ? new(lu, ipiv) : error("LU: dimension mismatch") + m == numel(ipiv) ? new(lu, ipiv) : error("LUDense: dimension mismatch") end end +size(A::LUDense) = size(A.lu) +size(A::LUDense,n) = size(A.lu,n) -function LU{T<:Union(Float64,Float32,Complex64,Complex128)}(A::Matrix{T}) - lu, ipiv = _jl_lapack_getrf(copy(A)) - LU{T}(lu, ipiv) +lu{T<:LapackScalar}(A::Matrix{T}) = lu!(copy(A)) +function lu!{T<:LapackScalar}(A::Matrix{T}) + lu, ipiv = _jl_lapack_getrf(A) + LUDense{T}(lu, ipiv) end -LU{T<:Number}(A::Matrix{T}) = LU(float64(A)) +lu{T<:Real}(A::Matrix{T}) = lu(float64(A)) -function det(lu::LU) +function det(lu::LUDense) m, n = size(lu.lu) if m != n error("det only defined for square matrices") end prod(diag(lu.lu)) * (bool(sum(lu.ipiv .!= 1:n) % 2) ? -1 : 1) end -det(A::Matrix) = det(LU(A)) +det(A::Matrix) = det(lu(A)) -(\){T<:Union(Float64,Float32,Complex128,Complex64)}(lu::LU{T}, B::StridedVecOrMat{T}) = +(\){T<:LapackScalar}(lu::LUDense{T}, B::StridedVecOrMat{T}) = _jl_lapack_getrs('N', lu.lu, lu.ipiv, copy(B)) -inv(lu::LU) = _jl_lapack_getri(copy(lu.lu), lu.ipiv) + +inv{T<:LapackScalar}(lu::LUDense{T}) = _jl_lapack_getri(copy(lu.lu), lu.ipiv) + +function factors{T<:LapackScalar}(lu::LUDense{T}) + LU, ipiv = lu.lu, lu.ipiv + m, n = size(LU) + + L = m >= n ? tril(LU, -1) + eye(m,n) : tril(LU, -1)[:, 1:m] + eye(m,m) + U = m <= n ? triu(LU) : triu(LU)[1:n, :] + P = [1:m] + for i=1:min(m,n) + t = P[i] + P[i] = P[ipiv[i]] + P[ipiv[i]] = t + end + L, U, P +end ## Multiplication by Q or Q' from a QR factorization for (orm2r, elty) in @@ -193,74 +238,72 @@ for (orm2r, elty) in end end -type QR{T} <: Factorization{T} +abstract QRFactorization{T} <: Factorization{T} + +## QR decomposition without column pivots +type QRDense{T} <: QRFactorization{T} hh::Matrix{T} # Householder transformations and R tau::Vector{T} # Scalar factors of transformations - function QR(hh::Matrix{T}, tt::Vector{T}) + function QRDense(hh::Matrix{T}, tt::Vector{T}) numel(tt) == min(size(hh)) ? new(hh, tt) : error("QR: mismatched dimensions") end end +size(A::QRFactorization) = size(A.hh) +size(A::QRFactorization,n) = size(A.hh,n) -function QR(A::Matrix) +function qr{T<:LapackScalar}(A::StridedMatrix{T}) hh, tt = _jl_lapack_geqrf(copy(A)) - QR{typeof(A[1])}(hh, tt) + QRDense{eltype(A)}(hh, tt) end -size{T<:Union(Float64,Float32,Complex128,Complex64)}(A::QR{T}) = size(A.hh) -size{T<:Union(Float64,Float32,Complex128,Complex64)}(A::QR{T},n) = size(A.hh,n) +qr{T<:Real}(x::StridedMatrix{T}) = qr(float64(x)) + +function factors{T<:LapackScalar}(qrd::QRDense{T}) + aa, tau = qrd.hh, qrd.tau + R = triu(aa[1:min(size(qrd)),:]) + _jl_lapack_orgqr(aa, tau), R +end ## Multiplication by Q from the QR decomposition -function (*){T<:Union(Float64,Float32,Complex128,Complex64)}(A::QR{T}, - B::StridedVecOrMat{T}) +(*){T<:LapackScalar}(A::QRFactorization{T}, B::StridedVecOrMat{T}) = _jl_lapack_orm2r('L', 'N', A.hh, size(A, 2), A.tau, copy(B)) -end ## Multiplication by Q' from the QR decomposition -function Ac_mul_B{T<:Union(Float64,Float32,Complex128,Complex64)}(A::QR{T}, - B::StridedVecOrMat{T}) +Ac_mul_B{T<:LapackScalar}(A::QRFactorization{T}, B::StridedVecOrMat{T}) = _jl_lapack_orm2r('L', iscomplex(A.tau) ? 'C' : 'T', A.hh, size(A, 2), A.tau, copy(B)) -end ## Least squares solution. Should be more careful about cases with m < n -function (\){T<:Union(Float64,Float32,Complex128,Complex64)}(A::QR{T}, - B::StridedVecOrMat{T}) +function (\){T<:LapackScalar}(A::QRDense{T}, B::StridedVecOrMat{T}) n = numel(A.tau) qtb = isa(B, Vector) ? (A' * B)[1:n] : (A' * B)[1:n, :] ## Not sure if this avoids copying A.hh[1:n,:] but at least it is not all of A.hh _jl_lapack_trtrs('U','N','N', A.hh[1:n,:], qtb) end -type QRP{T} <: Factorization{T} +type QRPDense{T} <: QRFactorization{T} hh::Matrix{T} tau::Vector{T} jpvt::Vector{Int32} - function QRP(hh::Matrix{T}, tt::Vector{T}, jj::Vector{Int32}) + function QRPDense(hh::Matrix{T}, tt::Vector{T}, jj::Vector{Int32}) m, n = size(hh) - numel(tt) == min(m,n) && numel(jj) == n ? new(hh,tt,jj) : error("QRP: mismatched dimensions") + numel(tt) == min(m,n) && numel(jj) == n ? new(hh,tt,jj) : error("QRPDense: mismatched dimensions") end end -function QRP(A::Matrix) - hh, tt, jj = _jl_lapack_geqp3(copy(A)) - QRP{typeof(A[1])}(hh, tt, jj) -end - -## Not sure how to avoid cut-and-paste programming on these. -## Create another abstract type with both QR{T} and QRP{T} as subtypes? -size{T<:Union(Float64,Float32,Complex128,Complex64)}(A::QRP{T}) = size(A.hh) -size{T<:Union(Float64,Float32,Complex128,Complex64)}(A::QRP{T},n) = size(A.hh,n) - -function (*){T<:Union(Float64,Float32,Complex128,Complex64)}(A::QRP{T}, - B::StridedVecOrMat{T}) - _jl_lapack_orm2r('L', 'N', A.hh, size(A, 2), A.tau, copy(B)) +function qrp{T<:LapackScalar}(A::StridedMatrix{T}) + aa, tau, jpvt = _jl_lapack_geqp3(copy(A)) + QRPDense{T}(aa, tau, jpvt) end +qrp{T<:Real}(x::StridedMatrix{T}) = qrp(float64(x)) -function Ac_mul_B{T<:Union(Float64,Float32,Complex128,Complex64)}(A::QRP{T}, - B::StridedVecOrMat{T}) - _jl_lapack_orm2r('L', iscomplex(A.tau) ? 'C' : 'T', A.hh, size(A, 2), A.tau, copy(B)) +function factors{T<:LapackScalar}(qrpd::QRPDense{T}) + aa, tau = qrpd.hh, qrpd.tau + R = triu(aa[1:min(size(qrpd)),:]) + _jl_lapack_orgqr(aa, tau), R, qrpd.jpvt end -function (\){T<:Union(Float64,Float32,Complex128,Complex64)}(A::QRP{T}, +# FIXME \ +function (\){T<:Union(Float64,Float32,Complex128,Complex64)}(A::QRPDense{T}, B::StridedVecOrMat{T}) n = numel(A.tau) ## Replace this with a direct call to _jl_lapack_trtrs to save copying A.hh? @@ -268,7 +311,7 @@ function (\){T<:Union(Float64,Float32,Complex128,Complex64)}(A::QRP{T}, triu(A.hh[1:n,:]) \ (A' * B)[1:n] end -function (\){T<:Union(Float64,Float32,Complex128,Complex64)}(A::QRP{T}, +function (\){T<:Union(Float64,Float32,Complex128,Complex64)}(A::QRPDense{T}, B::StridedVecOrMat{T}) ## may be better to define one method for B::Vector{T} and another for StridedMatrix BV = isa(B, Vector) @@ -288,13 +331,12 @@ type LDLTTridiagonal{T} <: Factorization{T} D::Vector{T} E::Vector{T} end -function LDLTTridiagonal{T<:LapackScalar}(A::Tridiagonal{T}) +function ldlt{T<:LapackScalar}(A::Tridiagonal{T}) D = copy(A.d) E = copy(A.dl) _jl_lapack_pttrf(D, E) - LDLTTridiagonal(D, E) + LDLTTridiagonal{T}(D, E) end -LDLT(A::Tridiagonal) = LDLTTridiagonal(A) (\){T<:LapackScalar}(C::LDLTTridiagonal{T}, B::StridedVecOrMat{T}) = _jl_lapack_pttrs(C.D, C.E, copy(B)) @@ -304,25 +346,21 @@ type LUTridiagonal{T} <: Factorization{T} ipiv::Vector{Int32} function LUTridiagonal(lu::Tridiagonal{T}, ipiv::Vector{Int32}) m, n = size(lu) - m == numel(ipiv) ? new(lu, ipiv) : error("LU: dimension mismatch") + m == numel(ipiv) ? new(lu, ipiv) : error("LUTridiagonal: dimension mismatch") end end show(io, lu::LUTridiagonal) = print(io, "LU decomposition of ", summary(lu.lu)) -function LU{T<:LapackScalar}(A::Tridiagonal{T}) +function lu{T<:LapackScalar}(A::Tridiagonal{T}) lu, ipiv = _jl_lapack_gttrf(copy(A)) LUTridiagonal{T}(lu, ipiv) end -function lu(A::Tridiagonal) - error("lu(A) is not defined when A is Tridiagonal. Use LU(A) instead.") -end - function det(lu::LUTridiagonal) prod(lu.lu.d) * (bool(sum(lu.ipiv .!= 1:n) % 2) ? -1 : 1) end -det(A::Tridiagonal) = det(LU(A)) +det(A::Tridiagonal) = det(lu(A)) (\){T<:LapackScalar}(lu::LUTridiagonal{T}, B::StridedVecOrMat{T}) = _jl_lapack_gttrs('N', lu.lu, lu.ipiv, copy(B)) diff --git a/base/linalg.jl b/base/linalg.jl index c68fefd0f799d..905e25e7a0bc0 100644 --- a/base/linalg.jl +++ b/base/linalg.jl @@ -23,6 +23,8 @@ triu(M::AbstractMatrix) = triu(M,0) tril(M::AbstractMatrix) = tril(M,0) #triu{T}(M::AbstractMatrix{T}, k::Integer) #tril{T}(M::AbstractMatrix{T}, k::Integer) +triu!(M::AbstractMatrix) = triu!(M,0) +tril!(M::AbstractMatrix) = tril!(M,0) #diff(a::AbstractVector) #diff(a::AbstractMatrix, dim::Integer) diff --git a/base/linalg_dense.jl b/base/linalg_dense.jl index af121e2940909..57ed7433a3873 100644 --- a/base/linalg_dense.jl +++ b/base/linalg_dense.jl @@ -264,6 +264,26 @@ triu{T}(M::Matrix{T}, k::Integer) = [ j-i >= k ? M[i,j] : zero(T) for i=1:size(M,1), j=1:size(M,2) ] tril{T}(M::Matrix{T}, k::Integer) = [ j-i <= k ? M[i,j] : zero(T) for i=1:size(M,1), j=1:size(M,2) ] +function triu!{T}(M::Matrix{T}, k::Integer) + m, n = size(M) + for i = 1:m + for j = 1:n + if j-i < k + M[i,j] = zero(T) + end + end + end +end +function tril!{T}(M::Matrix{T}, k::Integer) + m, n = size(M) + for i = 1:m + for j = 1:n + if j-i > k + M[i,j] = zero(T) + end + end + end +end diff(a::Vector) = [ a[i+1] - a[i] for i=1:length(a)-1 ] diff --git a/base/linalg_lapack.jl b/base/linalg_lapack.jl index db655ef84dace..ce5bb1960aa77 100644 --- a/base/linalg_lapack.jl +++ b/base/linalg_lapack.jl @@ -201,67 +201,6 @@ for (orgqr, elty) in end end -# chol() does not check that input matrix is symmetric/hermitian -# It simply uses upper triangular half -chol{T<:Integer}(x::StridedMatrix{T}) = chol(float64(x)) - -function chol{T<:Union(Float32,Float64,Complex64,Complex128)}(A::StridedMatrix{T}) - R = chol!(copy(A)) -end - -## chol! overwrites A with either the upper or lower triangular -## Cholesky factor (default is upper) -chol!{T<:Union(Float32,Float64,Complex64,Complex128)}(A::StridedMatrix{T}) = chol!(A, 'U') -function chol!{T<:Union(Float32,Float64,Complex64,Complex128)}(A::StridedMatrix{T}, uplo::LapackChar) - info = _jl_lapack_potrf(uplo, A) - if info != 0 error("chol: matrix is not positive definite, error $info") end - uplo == 'U' ? triu(A) : tril(A) -end - - -lu{T<:Integer}(x::StridedMatrix{T}) = lu(float64(x)) - -## LU decomposition returning L and U separately and P as a permutation -function lu{T<:Union(Float32,Float64,Complex64,Complex128)}(A::StridedMatrix{T}) - LU, ipiv = _jl_lapack_getrf(copy(A)) - m, n = size(A) - - L = m >= n ? tril(LU, -1) + eye(m,n) : tril(LU, -1)[:, 1:m] + eye(m,m) - U = m <= n ? triu(LU) : triu(LU)[1:n, :] - P = [1:m] - for i=1:min(m,n) - t = P[i] - P[i] = P[ipiv[i]] - P[ipiv[i]] = t - end - L, U, P -end - -## lu! overwrites A with the components of the LU decomposition -## Note that returned value includes the pivot indices, not the permutation -function lu!{T<:Union(Float32,Float64,Complex64,Complex128)}(A::StridedMatrix{T}) - _jl_lapack_getrf(A) -end - -## QR decomposition without column pivots -qr{T<:Integer}(x::StridedMatrix{T}) = qr(float64(x)) - -function qr{T<:Union(Float32,Float64,Complex64,Complex128)}(A::StridedMatrix{T}) - aa, tau = _jl_lapack_geqrf(copy(A)) - R = triu(aa[1:min(size(A)),:]) - _jl_lapack_orgqr(aa, tau), R -end - -## QR decomposition with column pivots -qrp{T<:Integer}(x::StridedMatrix{T}) = qrp(float64(x)) - -function qrp{T<:Union(Float32,Float64,Complex64,Complex128)}(A::StridedMatrix{T}) - aa, tau, jpvt = _jl_lapack_geqp3(copy(A)) - R = triu(aa[1:min(size(A)),:]) - _jl_lapack_orgqr(aa, tau), R, jpvt -end - - # eigenvalue-eigenvector, symmetric (Hermitian) or general cases for (syev, geev, elty) in ((:dsyev_,:dgeev_,:Float64), diff --git a/test/factorizations.jl b/test/factorizations.jl index 1b78ecdcfd9a6..3d17621007564 100644 --- a/test/factorizations.jl +++ b/test/factorizations.jl @@ -15,15 +15,15 @@ begin bd = rand(n) bz = complex(bd) # Cholesky decomposition - chd = Cholesky(symd) - chz = Cholesky(herz) + chd = chol(symd) + chz = chol(herz) xd = chd \ bd xz = chz \ bz @assert norm(symd * xd - bd) < Eps @assert norm(herz * xz - bz) < Eps # LU decomposition - lud = LU(sqd) - luz = LU(sqz) + lud = lu(sqd) + luz = lu(sqz) xd = lud \ bd xz = luz \ bz @assert norm(sqd * xd - bd) < Eps @@ -38,10 +38,10 @@ begin @assert_approx_eq det(symd) prod(eig(symd)[1]) # @assert_approx_eq det(herz) prod(eig(herz)[1]) - qrd = QR(mmd) # QR and QRP decompositions - qrz = QR(mmz) - qrpd = QRP(mmd) - qrpz = QRP(mmz) + qrd = qr(mmd) # QR and QRP decompositions + qrz = qr(mmz) + qrpd = qrp(mmd) + qrpz = qrp(mmz) yyd = randn(3*n) yyz = complex(yyd) qyd = qrd * yyd diff --git a/test/lapack.jl b/test/lapack.jl index 37c2d7a1e3b7e..b15777212c773 100644 --- a/test/lapack.jl +++ b/test/lapack.jl @@ -7,20 +7,20 @@ n = 10 a = rand(n,n) asym = a+a'+n*eye(n) b = rand(n) -r = chol(asym) # Cholesky decomposition +r = factors(chol(asym)) # Cholesky decomposition @assert norm(r'*r - asym) < Eps -l = chol!(copy(asym), 'L') # lower-triangular Cholesky decomposition +l = factors(chol!(copy(asym), 'L')) # lower-triangular Cholesky decomposition @assert norm(l*l' - asym) < Eps -(l,u,p) = lu(a) # LU decomposition +(l,u,p) = factors(lu(a)) # LU decomposition @assert norm(l*u - a[p,:]) < Eps @assert norm(l[invperm(p),:]*u - a) < Eps -(q,r) = qr(a) # QR decomposition +(q,r) = factors(qr(a)) # QR decomposition @assert norm(q*r - a) < Eps -(q,r,p) = qrp(a) # pivoted QR decomposition +(q,r,p) = factors(qrp(a)) # pivoted QR decomposition @assert norm(q*r - a[:,p]) < Eps @assert norm(q*r[:,invperm(p)] - a) < Eps @@ -166,7 +166,7 @@ invFv = F\v @assert norm(solve(T,v) - invFv) < Eps B = randn(n,2) @assert norm(solve(T, B) - F\B) < Eps -Tlu = LU(T) +Tlu = lu(T) x = Tlu\v @assert norm(x - invFv) < Eps @@ -174,7 +174,7 @@ x = Tlu\v Ts = Tridiagonal(dl, d, dl) Fs = full(Ts) invFsv = Fs\v -Tldlt = LDLT(Ts) +Tldlt = ldlt(Ts) x = Tldlt\v @assert norm(x - invFsv) < Eps