diff --git a/base/linalg/bidiag.jl b/base/linalg/bidiag.jl index 3190ca4032078..9c447f1b11e92 100644 --- a/base/linalg/bidiag.jl +++ b/base/linalg/bidiag.jl @@ -117,36 +117,43 @@ end #Linear solvers -\{T<:Number}(A::Union(Bidiagonal{T},Triangular{T}), b::AbstractVector{T}) = naivesub!(A, b, similar(b)) -\{T<:Number}(A::Union(Bidiagonal{T},Triangular{T}), B::AbstractMatrix{T}) = hcat([naivesub!(A, B[:,i], similar(B[:,i])) for i=1:size(B,2)]...) +function \{TA<:Number,Tb<:Number}(A::Union(Bidiagonal{TA},Triangular{TA}), b::AbstractVector{Tb}) + TAb = typeof(one(TA)/one(Tb)) + naivesub!(A, b, similar(b, TAb)) +end +function \{TA<:Number,TB<:Number}(A::Union(Bidiagonal{TA},Triangular{TA}), B::AbstractMatrix{TB}) + TAB = typeof(one(TA)/one(TB)) + hcat([naivesub!(A, B[:,i], similar(B[:,i], TAB)) for i=1:size(B,2)]...) +end A_ldiv_B!(A::Union(Bidiagonal,Triangular), b::AbstractVector)=naivesub!(A,b) At_ldiv_B!(A::Union(Bidiagonal,Triangular), b::AbstractVector)=naivesub!(transpose(A),b) Ac_ldiv_B!(A::Union(Bidiagonal,Triangular), b::AbstractVector)=naivesub!(ctranspose(A),b) for func in (:A_ldiv_B!, :Ac_ldiv_B!, :At_ldiv_B!) @eval begin - function ($func)(A::Bidiagonal, B::AbstractMatrix) + function ($func)(A::Union(Bidiagonal,Triangular), B::AbstractMatrix) + tmp = similar(B[:,1]) + n = size(B, 1) for i=1:size(B,2) - ($func)(A, B[:,i]) + copy!(tmp, 1, B, (i-1)*n+1, n) + ($func)(A, tmp) + copy!(B, (i-1)*n+1, tmp, 1, n) # Modify this when array view are implemented. end B end end end for func in (:A_ldiv_Bt!, :Ac_ldiv_Bt!, :At_ldiv_Bt!) @eval begin - function ($func)(A::Bidiagonal, B::AbstractMatrix) + function ($func)(A::Union(Bidiagonal,Triangular), B::AbstractMatrix) + tmp = similar(B[:,2]) + m, n = size(B) + nm = n*m for i=1:size(B,1) - ($func)(A, B[i,:]) + copy!(tmp, 1, B, i:m:nm, n) + ($func)(A, tmp) + copy!(B, i:m:nm, tmp, 1, n) end B end end end -function inv(A::Union(Bidiagonal, Triangular)) - B = eye(eltype(A), size(A, 1)) - for i=1:size(B,2) - naivesub!(A, B[:,i]) - end - B -end - #Generic solver using naive substitution function naivesub!{T}(A::Bidiagonal{T}, b::AbstractVector, x::AbstractVector=b) N = size(A, 2) diff --git a/base/linalg/dense.jl b/base/linalg/dense.jl index c2863a3e1c6dc..71af99b5412c1 100644 --- a/base/linalg/dense.jl +++ b/base/linalg/dense.jl @@ -412,12 +412,24 @@ function factorize!{T}(A::Matrix{T}) return Triangular(A, :U) end if herm - C, info = T <: BlasFloat ? LAPACK.potrf!('U', copy(A)) : LAPACK.potrf!('U', float(A)) + if T <: BlasFloat + C, info = LAPACK.potrf!('U', copy(A)) + elseif typeof(one(T)/one(T)) <: BlasFloat + C, info = LAPACK.potrf!('U', float(A)) + else + error("Unable to factorize hermitian $(typeof(A)). Try converting to other element type or use explicit factorization.") + end if info == 0 return Cholesky(C, 'U') end return factorize!(Hermitian(A)) end if sym - C, info = T <: BlasFloat ? LAPACK.potrf!('U', copy(A)) : LAPACK.potrf!('U', float(A)) + if T <: BlasFloat + C, info = LAPACK.potrf!('U', copy(A)) + elseif eltype(one(T)/one(T)) <: BlasFloat + C, info = LAPACK.potrf!('U', float(A)) + else + error("Unable to factorize symmetric $(typeof(A)). Try converting to other element type or use explicit factorization.") + end if info == 0 return Cholesky(C, 'U') end return factorize!(Symmetric(A)) end @@ -428,13 +440,8 @@ end factorize(A::AbstractMatrix) = factorize!(copy(A)) -(\){T1<:BlasFloat, T2<:BlasFloat}(A::StridedMatrix{T1}, B::StridedVecOrMat{T2}) = - (\)(convert(Array{promote_type(T1,T2)},A), convert(Array{promote_type(T1,T2)},B)) -(\){T1<:BlasFloat, T2<:Number}(A::StridedMatrix{T1}, B::StridedVecOrMat{T2}) = (\)(A, convert(Array{T1}, B)) -(\){T1<:Number, T2<:BlasFloat}(A::StridedMatrix{T1}, B::StridedVecOrMat{T2}) = (\)(convert(Array{T2}, A), B) -(\){T1<:Number, T2<:Number}(A::StridedMatrix{T1}, B::StridedVecOrMat{T2}) = T1<:Complex||T2<:Complex ? (\)(complex128(A), complex128(B)) : (\)(float64(A), float64(B)) (\)(a::Vector, B::StridedVecOrMat) = (\)(reshape(a, length(a), 1), B) -function (\){T<:BlasFloat}(A::StridedMatrix{T}, B::StridedVecOrMat{T}) +function (\)(A::StridedMatrix, B::StridedVecOrMat) m, n = size(A) if m == n if istril(A) diff --git a/base/linalg/diagonal.jl b/base/linalg/diagonal.jl index c75fd98334e70..8ccd36a871605 100644 --- a/base/linalg/diagonal.jl +++ b/base/linalg/diagonal.jl @@ -5,6 +5,8 @@ type Diagonal{T} <: AbstractMatrix{T} end Diagonal(A::Matrix) = Diagonal(diag(A)) +convert{T<:Number,S<:Number}(::Type{Diagonal{T}}, A::Diagonal{S}) = T == S ? A : Diagonal(convert(Vector{T}, A.diag)) + size(D::Diagonal) = (length(D.diag),length(D.diag)) size(D::Diagonal,d::Integer) = d<1 ? error("dimension out of range") : (d<=2 ? length(D.diag) : 1) diff --git a/base/linalg/factorization.jl b/base/linalg/factorization.jl index 786190d869c0e..b73a6c3bdbbfe 100644 --- a/base/linalg/factorization.jl +++ b/base/linalg/factorization.jl @@ -10,13 +10,6 @@ macro assertnonsingular(A, info) :(($info)==0 ? $A : throw(SingularException($info))) 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)) -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 # ########################## @@ -51,8 +44,8 @@ cholfact(x::Number) = @assertposdef Cholesky(fill(sqrt(x), 1, 1), :U) !(imag(x) chol(A::Union(Number, AbstractMatrix), uplo::Symbol) = cholfact(A, uplo)[uplo] chol(A::Union(Number, AbstractMatrix)) = triu!(cholfact(A, :U).UL) -size(C::Cholesky) = size(C.UL) -size(C::Cholesky,d::Integer) = size(C.UL,d) +size(C::Union(Cholesky, CholeskyPivoted)) = size(C.UL) +size(C::Union(Cholesky, CholeskyPivoted), d::Integer) = size(C.UL,d) function getindex(C::Cholesky, d::Symbol) d == :U && return triu!(symbol(C.uplo) == d ? C.UL : C.UL') @@ -126,24 +119,42 @@ rank(C::CholeskyPivoted) = C.rank #################### # LU Factorization # #################### -type LU{T<:BlasFloat} <: Factorization{T} +immutable LU{T} <: Factorization{T} factors::Matrix{T} ipiv::Vector{BlasInt} info::BlasInt end -LU{T<:BlasFloat}(A::StridedMatrix{T}) = LU{T}(LAPACK.getrf!(A)...) -lufact!(A::StridedMatrix) = lufact!(float(A)) -lufact!{T<:BlasFloat}(A::StridedMatrix{T}) = LU(A) -lufact{T<:BlasFloat}(A::StridedMatrix{T}) = lufact!(copy(A)) -lufact(A::StridedMatrix) = lufact!(float(A)) -lufact(x::Number) = LU(fill(x, 1, 1), [1], x == 0 ? 1 : 0) +lufact!{T<:BlasFloat}(A::StridedMatrix{T}) = LU(LAPACK.getrf!(A)...) +function lufact!{T}(A::AbstractMatrix{T}) + typeof(one(T)/one(T)) <: BlasFloat && return lufact!(float(A)) + m, n = size(A) + minmn = min(m,n) + info = 0 + for k = 1:minmn-1 + if A[k,k] == 0; info = k; break; end + for i = k+1:m + A[i,k] /= A[k,k] + end + for j = k+1:n + for i = k+1:m + A[i,j] -= A[i,k]*A[k,j] + end + end + end + if minmn > 0 && A[minmn,minmn] == 0; info = minmn; end + LU(A, BlasInt[1:minmn], convert(BlasInt, info)) +end +lufact(A::StridedMatrix) = lufact!(copy(A)) +lufact(x::Number) = LU(fill(x, 1, 1), BlasInt[1], x == 0 ? one(BlasInt) : zero(BlasInt)) function lu(A::Union(Number, AbstractMatrix)) F = lufact(A) F[:L], F[:U], F[:p] end +convert{T}(::Type{LU{T}}, F::LU) = LU(convert(Matrix{T}, F.factors), F.ipiv, F.info) + size(A::LU) = size(A.factors) size(A::LU,n) = size(A.factors,n) @@ -196,6 +207,7 @@ function logdet{T<:Complex}(A::LU{T}) end A_ldiv_B!{T<:BlasFloat}(A::LU{T}, B::StridedVecOrMat{T}) = @assertnonsingular LAPACK.getrs!('N', A.factors, A.ipiv, B) A.info +A_ldiv_B!(A::LU, B::StridedVecOrMat) = A_ldiv_B!(Triangular(A.factors, :U, false), A_ldiv_B!(Triangular(A.factors, :L, true), B)) 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 @@ -203,7 +215,7 @@ Ac_ldiv_Bc{T<:BlasComplex}(A::LU{T}, B::StridedVecOrMat{T}) = @assertnonsingular /{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 +inv{T<:BlasFloat}(A::LU{T})=@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))) @@ -213,13 +225,13 @@ cond(A::LU, p) = inv(LAPACK.gecon!(p == 1 ? '1' : 'I', A.factors, norm(A[:L][A[: # 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} + factors::Matrix{S} T::Matrix{S} end QR{T<:BlasFloat}(A::StridedMatrix{T}, nb::Integer = min(minimum(size(A)), 36)) = QR(LAPACK.geqrt!(A, nb)...) type QRPivoted{T} <: Factorization{T} - hh::Matrix{T} + factors::Matrix{T} tau::Vector{T} jpvt::Vector{BlasInt} end @@ -236,48 +248,57 @@ function qr(A::Union(Number, AbstractMatrix); pivot=false, thin::Bool=true) full(F[:Q], thin=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 == :R && return triu(A.factors[1:minimum(size(A)),:]) d == :Q && return QRPackedQ(A) throw(KeyError(d)) end type QRPackedQ{S} <: AbstractMatrix{S} - vs::Matrix{S} + factors::Matrix{S} T::Matrix{S} end -QRPackedQ(A::QR) = QRPackedQ(A.vs, A.T) +QRPackedQ(A::QR) = QRPackedQ(A.factors, A.T) +type QRPivotedQ{T} <: AbstractMatrix{T} + factors::Matrix{T} # Householder transformations and R + tau::Vector{T} # Scalar factors of transformations +end +QRPivotedQ(A::QRPivoted) = QRPivotedQ(A.factors, A.tau) -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) +size(A::Union(QR, QRPivoted), dim::Integer) = size(A.factors, dim) +size(A::Union(QR, QRPivoted)) = size(A.factors) +size(A::Union(QRPackedQ, QRPivotedQ), dim::Integer) = 0 < dim ? (dim <= 2 ? size(A.factors, 1) : 1) : throw(BoundsError()) +size(A::Union(QRPackedQ, QRPivotedQ)) = 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.factors)...) : eye(T, size(A.factors,1))) +function full{T<:BlasFloat}(A::QRPivotedQ{T}; thin::Bool=true) + m, n = size(A.factors) + Ahhpad = thin ? copy(A.factors) : [A.factors zeros(T, m, max(0, m - n))] + LAPACK.orgqr!(Ahhpad, A.tau) +end print_matrix(io::IO, A::QRPackedQ, rows::Integer, cols::Integer) = print_matrix(io, full(A, thin=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 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) + Bpad = size(B, 1)==size(A.factors, 2) ? [B; zeros(T, size(A.factors, 1) - size(A.factors, 2), size(B, 2))] : copy(B) + LAPACK.gemqrt!('L', 'N', A.factors, 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)) +*{T<:BlasFloat}(A::StridedVecOrMat{T}, B::QRPackedQ{T}) = LAPACK.gemqrt!('R', 'N', B.factors, B.T, copy(A)) +Ac_mul_B{T<:BlasReal}(A::QRPackedQ{T}, B::StridedVecOrMat{T}) = LAPACK.gemqrt!('L','T',A.factors,A.T,copy(B)) +Ac_mul_B{T<:BlasComplex}(A::QRPackedQ{T}, B::StridedVecOrMat{T}) = LAPACK.gemqrt!('L','C',A.factors,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) + Apad = size(A, 2)==size(B.factors, 2) ? [A zeros(T, size(A, 1), size(B.factors, 1) - size(B.factors, 2))] : copy(A) + LAPACK.gemqrt!('R', iseltype(B.factors,Complex) ? 'C' : 'T', B.factors, 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),:] -size(A::QRPivoted, args::Integer...) = size(A.hh, args...) - function getindex{T<:BlasFloat}(A::QRPivoted{T}, d::Symbol) - d == :R && return triu(A.hh[1:minimum(size(A)),:]) + d == :R && return triu(A.factors[1:minimum(size(A)),:]) d == :Q && return QRPivotedQ(A) d == :p && return A.jpvt if d == :P @@ -294,26 +315,26 @@ end # Julia implementation similarly to xgelsy function \{T<:BlasFloat}(A::QRPivoted{T}, B::StridedMatrix{T}, rcond::Real) - nr = minimum(size(A.hh)) + nr = minimum(size(A.factors)) nrhs = size(B, 2) if nr == 0 return zeros(0, nrhs), 0 end - ar = abs(A.hh[1]) + ar = abs(A.factors[1]) if ar == 0 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]) + tmin, smin, cmin = LAPACK.laic1!(2, sub(xmin, 1:rnk), tmin, sub(A.factors, 1:rnk, rnk + 1), A.factors[rnk + 1, rnk + 1]) + tmax, smax, cmax = LAPACK.laic1!(1, sub(xmax, 1:rnk), tmax, sub(A.factors, 1:rnk, rnk + 1), A.factors[rnk + 1, rnk + 1]) tmax*rcond > tmin && break xmin[1:rnk + 1] = [smin*sub(xmin, 1:rnk), cmin] xmax[1:rnk + 1] = [smax*sub(xmin, 1:rnk), cmax] 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)] + C, tau = LAPACK.tzrzf!(A.factors[1:rnk,:]) + X = [Triangular(C[1:rnk,1:rnk],:U)\(A[:Q]'B)[1:rnk,:]; zeros(T, size(A.factors, 2) - rnk, nrhs)] LAPACK.ormrz!('L', iseltype(B, Complex) ? 'C' : 'T', C, tau, X) return X[invperm(A[:p]),:], rnk end @@ -321,35 +342,19 @@ 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, thin=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) + Bpad = size(A.factors, 2)==size(B, 1) ? [B; zeros(T, size(A.factors, 1) - size(A.factors, 2), size(B, 2))] : copy(B) + LAPACK.ormqr!('L', 'N', A.factors, 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{T<:BlasReal}(A::QRPivotedQ{T}, B::StridedVecOrMat{T}) = LAPACK.ormqr!('L','T',A.factors,A.tau,copy(B)) +Ac_mul_B{T<:BlasComplex}(A::QRPivotedQ{T}, B::StridedVecOrMat{T}) = LAPACK.ormqr!('L','C',A.factors,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)) +*(A::StridedVecOrMat, B::QRPivotedQ) = LAPACK.ormqr!('R', 'N', B.factors, 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) + Apad = size(A, 2)==size(B.factors, 2) ? [A zeros(T, size(A, 1), size(B.factors, 1) - size(B.factors, 2))] : copy(A) + LAPACK.ormqr!('R', iseltype(B.factors,Complex) ? 'C' : 'T', B.factors, B.tau, Apad) end ##TODO: Add methods for rank(A::QRP{T}) and adjust the (\) method accordingly @@ -358,7 +363,7 @@ end # FIXME! Should add balancing option through xgebal type Hessenberg{T} <: Factorization{T} - hh::Matrix{T} + factors::Matrix{T} tau::Vector{T} function Hessenberg(hh::Matrix{T}, tau::Vector{T}) chksquare(hh) @@ -377,19 +382,19 @@ type HessenbergQ{T} <: AbstractMatrix{T} hh::Matrix{T} tau::Vector{T} end -HessenbergQ(A::Hessenberg) = HessenbergQ(A.hh, A.tau) -size(A::HessenbergQ, args...) = size(A.hh, args...) +HessenbergQ(A::Hessenberg) = HessenbergQ(A.factors, A.tau) +size(A::HessenbergQ, args...) = size(A.factors, args...) getindex(A::HessenbergQ, i::Real) = getindex(full(A), i) getindex(A::HessenbergQ, i::AbstractArray) = getindex(full(A), i) getindex(A::HessenbergQ, args...) = getindex(full(A), args...) function getindex(A::Hessenberg, d::Symbol) d == :Q && return HessenbergQ(A) - d == :H && return triu(A.hh, -1) + d == :H && return triu(A.factors, -1) throw(KeyError(d)) end -full(A::HessenbergQ) = LAPACK.orghr!(1, size(A.hh, 1), copy(A.hh), A.tau) +full(A::HessenbergQ) = LAPACK.orghr!(1, size(A.factors, 1), copy(A.factors), A.tau) # Eigenvalues type Eigen{T,V} <: Factorization{T} @@ -694,3 +699,20 @@ function schur(A::AbstractMatrix, B::AbstractMatrix) SchurF = schurfact(A, B) SchurF[:S], SchurF[:T], SchurF[:Q], SchurF[:Z] end + +### General promotion rules +inv{T}(F::Factorization{T}) = A_ldiv_B!(F, eye(T, size(F,1))) +function \{TF<:Number,TB<:Number}(F::Factorization{TF}, B::AbstractVecOrMat{TB}) + TFB = typeof(one(TF)/one(TB)) + A_ldiv_B!(convert(typeof(F).name.primary{TFB}, F), TB == TFB ? copy(B) : convert(typeof(B).name.primary{TFB}, B)) +end + +function Ac_ldiv_B{TF<:Number,TB<:Number}(F::Factorization{TF}, B::AbstractVecOrMat{TB}) + TFB = typeof(one(TF)/one(TB)) + Ac_ldiv_B!(convert(typeof(F).name.primary{TFB}, F), TB == TFB ? copy(B) : convert(typeof(B).name.primary{TFB}, B)) +end + +function At_ldiv_B{TF<:Number,TB<:Number}(F::Factorization{TF}, B::AbstractVecOrMat{TB}) + TFB = typeof(one(TF)/one(TB)) + At_ldiv_B!(convert(typeof(F).name.primary{TFB}, F), TB == TFB ? copy(B) : convert(typeof(B).name.primary{TFB}, B)) +end diff --git a/base/linalg/generic.jl b/base/linalg/generic.jl index 775489ca9e2e3..11c22a481ca0d 100644 --- a/base/linalg/generic.jl +++ b/base/linalg/generic.jl @@ -65,22 +65,37 @@ end norm{T<:Integer}(x::AbstractVector{T}, p::Number) = norm(float(x), p) norm(x::AbstractVector) = norm(x, 2) -function norm(A::AbstractMatrix, p::Number=2) +function norm{T}(A::AbstractMatrix{T}, p::Number=2) m, n = size(A) if m == 0 || n == 0 - a = zero(eltype(A)) + return zero(real(zero(T))) elseif m == 1 || n == 1 - a = norm(reshape(A, length(A)), p) + return norm(reshape(A, length(A)), p) elseif p == 1 - a = maximum(sum(abs(A),1)) + nrm = zero(real(zero(T))) + for j = 1:n + nrmj = zero(real(zero(T))) + for i = 1:m + nrmj += abs(A[i,j]) + end + nrm = max(nrm,nrmj) + end + return nrm elseif p == 2 - a = maximum(svdvals(A)) + return svdvals(A)[1] elseif p == Inf - a = maximum(sum(abs(A),2)) + nrm = zero(real(zero(T))) + for i = 1:m + nrmi = zero(real(zero(T))) + for j = 1:n + nrmi += abs(A[i,j]) + end + nrm = max(nrm,nrmi) + end + return nrm else throw(ArgumentError("invalid p-norm p=$p. Valid: 1, 2, Inf")) end - float(a) end norm(x::Number, p=nothing) = abs(x) @@ -109,10 +124,11 @@ trace(x::Number) = x #det(a::AbstractMatrix) inv(a::AbstractVector) = error("argument must be a square matrix") +inv{T}(A::AbstractMatrix{T}) = A_ldiv_B!(A,eye(T, chksquare(A))) function \{TA<:Number,TB<:Number}(A::AbstractMatrix{TA}, B::AbstractVecOrMat{TB}) TC = typeof(one(TA)/one(TB)) - return TB == TC ? A_ldiv_B!(A, copy(B)) : A_ldiv_B!(A, convert(Array{TC}, B)) + A_ldiv_B!(convert(typeof(A).name.primary{TC}, A), TB == TC ? copy(B) : convert(typeof(B).name.primary{TC}, B)) end \(a::AbstractVector, b::AbstractArray) = reshape(a, length(a), 1) \ b /(A::AbstractVecOrMat, B::AbstractVecOrMat) = (B' \ A')' diff --git a/base/linalg/triangular.jl b/base/linalg/triangular.jl index bbbd3d7980ea5..73f7866019e71 100644 --- a/base/linalg/triangular.jl +++ b/base/linalg/triangular.jl @@ -1,5 +1,5 @@ ## Triangular -type Triangular{T<:Number} <: AbstractMatrix{T} +immutable Triangular{T<:Number} <: AbstractMatrix{T} UL::Matrix{T} uplo::Char unitdiag::Char @@ -89,8 +89,8 @@ getindex{T}(A::Triangular{T}, i::Integer, j::Integer) = i == j ? (A.unitdiag == istril(A::Triangular) = A.uplo == 'L' || istriu(A.UL) istriu(A::Triangular) = A.uplo == 'U' || istril(A.UL) -transpose(A::Triangular) = Triangular(A.UL, A.uplo=='U':'L':'U', A.unitdiag) -ctranspose(A::Triangular) = conj(transpose(A)) +transpose(A::Triangular) = Triangular(symmetrize!(A.UL, A.uplo), A.uplo=='U'?'L':'U', A.unitdiag) +ctranspose(A::Triangular) = Triangular(symmetrize_conj!(A.UL, A.uplo), A.uplo=='U'?'L':'U', A.unitdiag) diag(A::Triangular) = diag(A.UL) big(A::Triangular) = Triangular(big(A.UL), A.uplo, A.unitdiag) diff --git a/base/sparse/sparsematrix.jl b/base/sparse/sparsematrix.jl index 7eb1a06a2658c..a14ceb8f6d9c3 100644 --- a/base/sparse/sparsematrix.jl +++ b/base/sparse/sparsematrix.jl @@ -150,6 +150,17 @@ function convert{Tv,Ti,TvS,TiS}(::Type{SparseMatrixCSC{Tv,Ti}}, S::SparseMatrixC end end +function convert{Tv,TvS,TiS}(::Type{SparseMatrixCSC{Tv}}, S::SparseMatrixCSC{TvS,TiS}) + if Tv == TvS + return S + else + return SparseMatrixCSC(S.m, S.n, + S.colptr, + S.rowval, + convert(Vector{Tv},S.nzval)) + end +end + function convert{Tv,Ti}(::Type{SparseMatrixCSC{Tv,Ti}}, M::Matrix) m, n = size(M) (I, J, V) = findnz(M) diff --git a/base/test.jl b/base/test.jl index 0bfa7b1d0c574..666ff7f032594 100644 --- a/base/test.jl +++ b/base/test.jl @@ -73,7 +73,7 @@ function test_approx_eq(va, vb, Eps, astr, bstr) end test_approx_eq(va, vb, astr, bstr) = - test_approx_eq(va, vb, 10^4*length(va)*eps(float(max(maximum(abs(va)), maximum(abs(vb))))), astr, bstr) + test_approx_eq(va, vb, 10^4*length(va)*max(eps(float(maximum(abs(va)))), eps(float(maximum(abs(vb))))), astr, bstr) macro test_approx_eq_eps(a, b, c) :(test_approx_eq($(esc(a)), $(esc(b)), $(esc(c)), $(string(a)), $(string(b)))) diff --git a/test/bitarray.jl b/test/bitarray.jl index c0f94d1ad0344..5ad1dbb982ce9 100644 --- a/test/bitarray.jl +++ b/test/bitarray.jl @@ -288,8 +288,8 @@ end b2 = randbool(n1, n1) @check_bit_operation (*) Matrix{Int} (b1, b2) -@check_bit_operation (/) Matrix{Float64} (b1, b1) -@check_bit_operation (\) Matrix{Float64} (b1, b1) +@check_bit_operation (/) Matrix{Float32} (b1, b1) +@check_bit_operation (\) Matrix{Float32} (b1, b1) b0 = falses(0) @check_bit_operation (&) BitVector (b0, b0) diff --git a/test/linalg.jl b/test/linalg.jl index be60400354dde..16e02c6e0679b 100644 --- a/test/linalg.jl +++ b/test/linalg.jl @@ -16,8 +16,8 @@ end areal = randn(n,n) aimg = randn(n,n) -breal = randn(n) -bimg = randn(n) +breal = randn(n,2) +bimg = randn(n,2) for elty in (Float32, Float64, Complex64, Complex128, Int) if elty == Complex64 || elty == Complex128 a = complex(areal, aimg) @@ -26,12 +26,12 @@ for elty in (Float32, Float64, Complex64, Complex128, Int) a = areal b = breal end - if elty != Int + if elty <: BlasFloat a = convert(Matrix{elty}, a) - b = convert(Vector{elty}, b) + b = convert(Matrix{elty}, b) else - a = rand(1:100, n, n) - b = rand(1:100, n) + a = rand(1:10, n, n) + b = rand(1:10, n, 2) end asym = a' + a # symmetric indefinite apd = a'*a # symmetric positive-definite @@ -569,10 +569,10 @@ end n=12 for relty in (Float16, Float32, Float64, BigFloat), elty in (relty, Complex{relty}) A = convert(Matrix{elty}, randn(n, n)) - b = convert(Vector{elty}, randn(n)) + b = convert(Matrix{elty}, randn(n, 2)) if elty <: Complex A += im*convert(Matrix{elty}, randn(n, n)) - b += im*convert(Vector{elty}, randn(n)) + b += im*convert(Matrix{elty}, randn(n, 2)) end for M in (triu(A), tril(A)) @@ -583,7 +583,7 @@ for relty in (Float16, Float32, Float64, BigFloat), elty in (relty, Complex{relt tx = TM \ b @test norm(x-tx,Inf) <= 4*condM*max(eps()*norm(tx,Inf), eps(relty)*norm(x,Inf)) if elty <: BlasFloat #test naivesub! against LAPACK - tx = LinAlg.naivesub!(TM, b, similar(b)) + tx = [LinAlg.naivesub!(TM, b[:,1]) LinAlg.naivesub!(TM, b[:,2])] @test norm(x-tx,Inf) <= 4*condM*max(eps()*norm(tx,Inf), eps(relty)*norm(x,Inf)) end @@ -668,11 +668,11 @@ end for relty in (Float16, Float32, Float64, BigFloat), elty in (relty, Complex{relty}) dv = convert(Vector{elty}, randn(n)) ev = convert(Vector{elty}, randn(n-1)) - b = convert(Vector{elty}, randn(n)) + b = convert(Matrix{elty}, randn(n, 2)) if (elty <: Complex) dv += im*convert(Vector{elty}, randn(n)) ev += im*convert(Vector{elty}, randn(n-1)) - b += im*convert(Vector{elty}, randn(n)) + b += im*convert(Matrix{elty}, randn(n, 2)) end for isupper in (true, false) #Test upper and lower bidiagonal matrices T = Bidiagonal(dv, ev, isupper) @@ -753,8 +753,6 @@ for relty in (Float16, Float32, Float64, BigFloat), elty in (relty, Complex{relt end #Test interconversion between special matrix types -using Base.Test - N=12 A=Diagonal([1:N]*1.0) for newtype in [Diagonal, Bidiagonal, SymTridiagonal, Tridiagonal, Triangular, Matrix] @@ -829,6 +827,26 @@ for elty in (Int32, Int64, Float32, Float64, Complex64, Complex128) @test_approx_eq gradient(x) g end +# Test rational matrices +## Integrate in general tests when more linear algebra is implemented in julia +a = convert(Matrix{Rational{BigInt}}, rand(1:10//1,n,n))/n +b = rand(1:10,n,2) +lua = factorize(a) +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(n) +@test_approx_eq a*(lua\b) b +@test_approx_eq det(a) det(float64(float(a))) +## Hilbert Matrix (very ill conditioned) +## Testing Rational{BigInt} and BigFloat version +nHilbert = 50 +H = Rational{BigInt}[1//(i+j-1) for i = 1:nHilbert,j = 1:nHilbert] +Hinv = Rational{BigInt}[(-1)^(i+j)*(i+j-1)*binomial(nHilbert+i-1,nHilbert-j)*binomial(nHilbert+j-1,nHilbert-i)*binomial(i+j-2,i-1)^2 for i = big(1):nHilbert,j=big(1):nHilbert] +@test inv(H) == Hinv +with_bigfloat_precision(2^10) do + @test norm(float64(inv(float(H)) - float(Hinv))) < 1e-100 +end ## Issue related tests # issue 1447 let