diff --git a/base/linalg/bidiag.jl b/base/linalg/bidiag.jl index 8701c11b8bdf5..bf8fd97418af7 100644 --- a/base/linalg/bidiag.jl +++ b/base/linalg/bidiag.jl @@ -114,7 +114,7 @@ end /(A::Bidiagonal, B::Number) = Bidiagonal(A.dv/B, A.ev/B, A.isupper) ==(A::Bidiagonal, B::Bidiagonal) = (A.dv==B.dv) && (A.ev==B.ev) && (A.isupper==B.isupper) -SpecialMatrix = Union(Diagonal, Bidiagonal, SymTridiagonal, Tridiagonal, TriangularUnion) +SpecialMatrix = Union(Diagonal, Bidiagonal, SymTridiagonal, Tridiagonal, AbstractTriangular) *(A::SpecialMatrix, B::SpecialMatrix)=full(A)*full(B) #Generic multiplication @@ -126,11 +126,11 @@ end #Linear solvers -A_ldiv_B!(A::Union(Bidiagonal, TriangularUnion), b::AbstractVector) = naivesub!(A, b) -At_ldiv_B!(A::Union(Bidiagonal, TriangularUnion), b::AbstractVector) = naivesub!(transpose(A), b) -Ac_ldiv_B!(A::Union(Bidiagonal, TriangularUnion), b::AbstractVector) = naivesub!(ctranspose(A), b) +A_ldiv_B!(A::Union(Bidiagonal, AbstractTriangular), b::AbstractVector) = naivesub!(A, b) +At_ldiv_B!(A::Union(Bidiagonal, AbstractTriangular), b::AbstractVector) = naivesub!(transpose(A), b) +Ac_ldiv_B!(A::Union(Bidiagonal, AbstractTriangular), b::AbstractVector) = naivesub!(ctranspose(A), b) for func in (:A_ldiv_B!, :Ac_ldiv_B!, :At_ldiv_B!) @eval begin - function ($func)(A::Union(Bidiagonal, TriangularUnion), B::AbstractMatrix) + function ($func)(A::Union(Bidiagonal, AbstractTriangular), B::AbstractMatrix) tmp = similar(B[:,1]) n = size(B, 1) for i = 1:size(B,2) @@ -142,7 +142,7 @@ for func in (:A_ldiv_B!, :Ac_ldiv_B!, :At_ldiv_B!) @eval begin end end end for func in (:A_ldiv_Bt!, :Ac_ldiv_Bt!, :At_ldiv_Bt!) @eval begin - function ($func)(A::Union(Bidiagonal, TriangularUnion), B::AbstractMatrix) + function ($func)(A::Union(Bidiagonal, AbstractTriangular), B::AbstractMatrix) tmp = similar(B[:, 2]) m, n = size(B) nm = n*m diff --git a/base/linalg/factorization.jl b/base/linalg/factorization.jl index dbf80951b8a08..deb3981ac8655 100644 --- a/base/linalg/factorization.jl +++ b/base/linalg/factorization.jl @@ -274,7 +274,7 @@ function A_mul_Bc!{T}(A::AbstractMatrix{T},Q::QRPackedQ{T}) end A end -A_mul_Bc(A::TriangularUnion, B::Union(QRCompactWYQ,QRPackedQ)) = A_mul_Bc(full(A), B) +A_mul_Bc(A::AbstractTriangular, B::Union(QRCompactWYQ,QRPackedQ)) = A_mul_Bc(full(A), B) function A_mul_Bc{TA,TB}(A::AbstractArray{TA}, B::Union(QRCompactWYQ{TB},QRPackedQ{TB})) TAB = promote_type(TA,TB) A_mul_Bc!(size(A,2)==size(B.factors,1) ? (TA == TAB ? copy(A) : convert(AbstractMatrix{TAB}, A)) : diff --git a/base/linalg/generic.jl b/base/linalg/generic.jl index 861d101afe83b..b54f4d9db6411 100644 --- a/base/linalg/generic.jl +++ b/base/linalg/generic.jl @@ -1,9 +1,5 @@ ## linalg.jl: Some generic Linear Algebra definitions -# Fall back arithmetic -+(A::AbstractMatrix, B::AbstractMatrix) = full(A) + full(B) --(A::AbstractMatrix, B::AbstractMatrix) = full(A) - full(B) - scale(X::AbstractArray, s::Number) = scale!(copy(X), s) scale(s::Number, X::AbstractArray) = scale!(copy(X), s) diff --git a/base/linalg/sparse.jl b/base/linalg/sparse.jl index d507bbb8e6a53..08a72a0f9b81e 100644 --- a/base/linalg/sparse.jl +++ b/base/linalg/sparse.jl @@ -126,7 +126,7 @@ end *{TvA,TiA}(X::BitArray{2}, A::SparseMatrixCSC{TvA,TiA}) = invoke(*, (AbstractMatrix, SparseMatrixCSC), X, A) # TODO: Tridiagonal * Sparse should be implemented more efficiently *{TX,TvA,TiA}(X::Tridiagonal{TX}, A::SparseMatrixCSC{TvA,TiA}) = invoke(*, (Tridiagonal, AbstractMatrix), X, A) -*{TvA,TiA}(X::TriangularUnion, A::SparseMatrixCSC{TvA,TiA}) = full(X)*A +*{TvA,TiA}(X::AbstractTriangular, A::SparseMatrixCSC{TvA,TiA}) = full(X)*A function *{TX,TvA,TiA}(X::AbstractMatrix{TX}, A::SparseMatrixCSC{TvA,TiA}) mX, nX = size(X) nX == A.m || throw(DimensionMismatch()) diff --git a/base/linalg/special.jl b/base/linalg/special.jl index 0c0af842fe4c0..0961652c7f8f1 100644 --- a/base/linalg/special.jl +++ b/base/linalg/special.jl @@ -46,12 +46,12 @@ function convert(::Type{SymTridiagonal}, A::Tridiagonal) SymTridiagonal(A.d, A.dl) end -function convert(::Type{Diagonal}, A::TriangularUnion) +function convert(::Type{Diagonal}, A::AbstractTriangular) full(A) == diagm(diag(A)) || throw(ArgumentError("Matrix cannot be represented as Diagonal")) Diagonal(diag(A)) end -function convert(::Type{Bidiagonal}, A::TriangularUnion) +function convert(::Type{Bidiagonal}, A::AbstractTriangular) fA = full(A) if fA == diagm(diag(A)) + diagm(diag(fA, 1), 1) return Bidiagonal(diag(A), diag(fA,1), true) @@ -62,9 +62,9 @@ function convert(::Type{Bidiagonal}, A::TriangularUnion) end end -convert(::Type{SymTridiagonal}, A::TriangularUnion) = convert(SymTridiagonal, convert(Tridiagonal, A)) +convert(::Type{SymTridiagonal}, A::AbstractTriangular) = convert(SymTridiagonal, convert(Tridiagonal, A)) -function convert(::Type{Tridiagonal}, A::TriangularUnion) +function convert(::Type{Tridiagonal}, A::AbstractTriangular) fA = full(A) if fA == diagm(diag(A)) + diagm(diag(fA, 1), 1) + diagm(diag(fA, -1), -1) return Tridiagonal(diag(fA, -1), diag(A), diag(fA,1)) @@ -115,7 +115,7 @@ for op in (:+, :-) end end -A_mul_Bc!(A::TriangularUnion, B::QRCompactWYQ) = A_mul_Bc!(full!(A),B) -A_mul_Bc!(A::TriangularUnion, B::QRPackedQ) = A_mul_Bc!(full!(A),B) -A_mul_Bc(A::TriangularUnion, B::Union(QRCompactWYQ,QRPackedQ)) = A_mul_Bc(full(A), B) +A_mul_Bc!(A::AbstractTriangular, B::QRCompactWYQ) = A_mul_Bc!(full!(A),B) +A_mul_Bc!(A::AbstractTriangular, B::QRPackedQ) = A_mul_Bc!(full!(A),B) +A_mul_Bc(A::AbstractTriangular, B::Union(QRCompactWYQ,QRPackedQ)) = A_mul_Bc(full(A), B) diff --git a/base/linalg/triangular.jl b/base/linalg/triangular.jl index 2ec3bb7881da7..086dacc1da353 100644 --- a/base/linalg/triangular.jl +++ b/base/linalg/triangular.jl @@ -1,9 +1,11 @@ ## Triangular +abstract AbstractTriangular{T,S<:AbstractMatrix} <: AbstractMatrix{T} # could be renamed to Triangular when than name has been fully deprecated + # First loop through all methods that don't need special care for upper/lower and unit diagonal for t in (:LowerTriangular, :UnitLowerTriangular, :UpperTriangular, :UnitUpperTriangular) @eval begin - immutable $t{T,S<:AbstractMatrix} <: AbstractMatrix{T} + immutable $t{T,S<:AbstractMatrix} <: AbstractTriangular{T,S} data::S end function $t(A::AbstractMatrix) @@ -35,12 +37,9 @@ for t in (:LowerTriangular, :UnitLowerTriangular, :UpperTriangular, :UnitUpperTr end end -# some cases can be handled with a union -typealias TriangularUnion{T,S} Union(UpperTriangular{T,S}, UnitUpperTriangular{T,S}, LowerTriangular{T,S}, UnitLowerTriangular{T,S}) - -full{T,S}(A::TriangularUnion{T,S}) = convert(Matrix, A) +full(A::AbstractTriangular) = convert(Matrix, A) -fill!(A::TriangularUnion, x) = (fill!(A.data, x); A) +fill!(A::AbstractTriangular, x) = (fill!(A.data, x); A) # then handle all methods that requires specific handling of upper/lower and unit diagonal @@ -102,7 +101,7 @@ function full!{T,S}(A::UnitUpperTriangular{T,S}) B end -getindex(A::TriangularUnion, i::Integer) = ((m, n) = divrem(i - 1, size(A, 1)); A[n + 1, m + 1]) +getindex(A::AbstractTriangular, i::Integer) = ((m, n) = divrem(i - 1, size(A, 1)); A[n + 1, m + 1]) getindex{T,S}(A::UnitLowerTriangular{T,S}, i::Integer, j::Integer) = i == j ? one(T) : (i > j ? A.data[i,j] : zero(A.data[i,j])) getindex{T,S}(A::LowerTriangular{T,S}, i::Integer, j::Integer) = i >= j ? A.data[i,j] : zero(A.data[i,j]) getindex{T,S}(A::UnitUpperTriangular{T,S}, i::Integer, j::Integer) = i == j ? one(T) : (i < j ? A.data[i,j] : zero(A.data[i,j])) @@ -172,6 +171,7 @@ end +(A::UnitLowerTriangular, B::LowerTriangular) = LowerTriangular(tril(A.data, -1) + B.data + I) +(A::UnitUpperTriangular, B::UnitUpperTriangular) = UpperTriangular(triu(A.data, 1) + triu(B.data, 1) + 2I) +(A::UnitLowerTriangular, B::UnitLowerTriangular) = LowerTriangular(tril(A.data, -1) + tril(B.data, -1) + 2I) ++(A::AbstractTriangular, B::AbstractTriangular) = full(A) + full(B) -(A::UpperTriangular, B::UpperTriangular) = UpperTriangular(A.data - B.data) -(A::LowerTriangular, B::LowerTriangular) = LowerTriangular(A.data - B.data) @@ -181,14 +181,15 @@ end -(A::UnitLowerTriangular, B::LowerTriangular) = LowerTriangular(tril(A.data, -1) - B.data + I) -(A::UnitUpperTriangular, B::UnitUpperTriangular) = UpperTriangular(triu(A.data, 1) - triu(B.data, 1)) -(A::UnitLowerTriangular, B::UnitLowerTriangular) = LowerTriangular(tril(A.data, -1) - tril(B.data, -1)) +-(A::AbstractTriangular, B::AbstractTriangular) = full(A) - full(B) ###################### # BlasFloat routines # ###################### -A_mul_B!(A::Tridiagonal, B::TriangularUnion) = A*full!(B) -A_mul_B!(C::AbstractVecOrMat, A::TriangularUnion, B::AbstractVecOrMat) = A_mul_B!(A, copy!(C, B)) -A_mul_Bc!(C::AbstractVecOrMat, A::TriangularUnion, B::AbstractVecOrMat) = A_mul_Bc!(A, copy!(C, B)) +A_mul_B!(A::Tridiagonal, B::AbstractTriangular) = A*full!(B) +A_mul_B!(C::AbstractVecOrMat, A::AbstractTriangular, B::AbstractVecOrMat) = A_mul_B!(A, copy!(C, B)) +A_mul_Bc!(C::AbstractVecOrMat, A::AbstractTriangular, B::AbstractVecOrMat) = A_mul_Bc!(A, copy!(C, B)) for (t, uploc, isunitc) in ((:LowerTriangular, 'L', 'N'), (:UnitLowerTriangular, 'L', 'U'), @@ -238,8 +239,8 @@ for (t, uploc, isunitc) in ((:LowerTriangular, 'L', 'N'), end end -errorbounds{T<:Union(BigFloat, Complex{BigFloat}),S<:StridedMatrix}(A::TriangularUnion{T,S}, X::StridedVecOrMat{T}, B::StridedVecOrMat{T}) = error("not implemented yet! Please submit a pull request.") -function errorbounds{TA<:Number,S<:StridedMatrix,TX<:Number,TB<:Number}(A::TriangularUnion{TA,S}, X::StridedVecOrMat{TX}, B::StridedVecOrMat{TB}) +errorbounds{T<:Union(BigFloat, Complex{BigFloat}),S<:StridedMatrix}(A::AbstractTriangular{T,S}, X::StridedVecOrMat{T}, B::StridedVecOrMat{T}) = error("not implemented yet! Please submit a pull request.") +function errorbounds{TA<:Number,S<:StridedMatrix,TX<:Number,TB<:Number}(A::AbstractTriangular{TA,S}, X::StridedVecOrMat{TX}, B::StridedVecOrMat{TB}) TAXB = promote_type(TA, TB, TX, Float32) errorbounds(convert(AbstractMatrix{TAXB}, A), convert(AbstractArray{TAXB}, X), convert(AbstractArray{TAXB}, B)) end @@ -733,12 +734,12 @@ end for f in (:*, :Ac_mul_B, :At_mul_B, :A_mul_Bc, :A_mul_Bt, :Ac_mul_Bc, :At_mul_Bt, :\, :Ac_ldiv_B, :At_ldiv_B) @eval begin - ($f)(A::TriangularUnion, B::TriangularUnion) = ($f)(A, full(B)) + ($f)(A::AbstractTriangular, B::AbstractTriangular) = ($f)(A, full(B)) end end for f in (:A_mul_Bc, :A_mul_Bt, :Ac_mul_Bc, :At_mul_Bt, :/, :A_rdiv_Bc, :A_rdiv_Bt) @eval begin - ($f)(A::TriangularUnion, B::TriangularUnion) = ($f)(full(A), B) + ($f)(A::AbstractTriangular, B::AbstractTriangular) = ($f)(full(A), B) end end @@ -746,7 +747,7 @@ end ### Multiplication with triangle to the left and hence rhs cannot be transposed. for (f, g) in ((:*, :A_mul_B!), (:Ac_mul_B, :Ac_mul_B!), (:At_mul_B, :At_mul_B!)) @eval begin - function ($f){TA,TB}(A::TriangularUnion{TA}, B::StridedVecOrMat{TB}) + function ($f){TA,TB}(A::AbstractTriangular{TA}, B::StridedVecOrMat{TB}) TAB = typeof(zero(TA)*zero(TB) + zero(TA)*zero(TB)) ($g)(TA == TAB ? copy(A) : convert(AbstractArray{TAB}, A), TB == TAB ? copy(B) : convert(AbstractArray{TAB}, B)) end @@ -789,7 +790,7 @@ end ### Multiplication with triangle to the rigth and hence lhs cannot be transposed. for (f, g) in ((:*, :A_mul_B!), (:A_mul_Bc, :A_mul_Bc!), (:A_mul_Bt, :A_mul_Bt!)) @eval begin - function ($f){TA,TB}(A::StridedVecOrMat{TA}, B::TriangularUnion{TB}) + function ($f){TA,TB}(A::StridedVecOrMat{TA}, B::AbstractTriangular{TB}) TAB = typeof(zero(TA)*zero(TB) + zero(TA)*zero(TB)) ($g)(TA == TAB ? copy(A) : convert(AbstractArray{TAB}, A), TB == TAB ? copy(B) : convert(AbstractArray{TAB}, B)) end @@ -866,8 +867,8 @@ sqrtm(A::LowerTriangular) = sqrtm(A.').' sqrtm(A::UnitLowerTriangular) = sqrtm(A.').' #Generic eigensystems -eigvals(A::TriangularUnion) = diag(A) -function eigvecs{T}(A::TriangularUnion{T}) +eigvals(A::AbstractTriangular) = diag(A) +function eigvecs{T}(A::AbstractTriangular{T}) TT = promote_type(T, Float32) TT <: BlasFloat && return eigvecs(convert(AbstractMatrix{TT}, A)) error("type not handled yet. Please submit a pull request.") @@ -877,14 +878,14 @@ det{T}(A::UnitLowerTriangular{T}) = one(T)*one(T) det{T}(A::UpperTriangular{T}) = prod(diag(A.data)) det{T}(A::LowerTriangular{T}) = prod(diag(A.data)) -eigfact(A::TriangularUnion) = Eigen(eigvals(A), eigvecs(A)) +eigfact(A::AbstractTriangular) = Eigen(eigvals(A), eigvecs(A)) #Generic singular systems for func in (:svd, :svdfact, :svdfact!, :svdvals) @eval begin - ($func)(A::TriangularUnion) = ($func)(full(A)) + ($func)(A::AbstractTriangular) = ($func)(full(A)) end end -factorize(A::TriangularUnion) = A +factorize(A::AbstractTriangular) = A diff --git a/base/linalg/uniformscaling.jl b/base/linalg/uniformscaling.jl index 4bf9627540511..9815a93bf2f23 100644 --- a/base/linalg/uniformscaling.jl +++ b/base/linalg/uniformscaling.jl @@ -84,7 +84,7 @@ inv(J::UniformScaling) = UniformScaling(inv(J.λ)) /(J::UniformScaling, x::Number) = UniformScaling(J.λ/x) \(J1::UniformScaling, J2::UniformScaling) = J1.λ == 0 ? throw(SingularException(1)) : UniformScaling(J1.λ\J2.λ) -\{T<:Number}(A::Union(Bidiagonal{T},TriangularUnion{T}), J::UniformScaling) = inv(A)*J.λ +\{T<:Number}(A::Union(Bidiagonal{T},AbstractTriangular{T}), J::UniformScaling) = inv(A)*J.λ \(J::UniformScaling, A::AbstractVecOrMat) = J.λ == 0 ? throw(SingularException(1)) : J.λ\A \(A::AbstractMatrix, J::UniformScaling) = inv(A)*J.λ