Skip to content

Commit

Permalink
Remove + and - fall backs for AbstractMatrix. Use abstract type inste…
Browse files Browse the repository at this point in the history
…ad of union.
  • Loading branch information
andreasnoack committed Jan 19, 2015
1 parent 56ffa13 commit e40c11b
Show file tree
Hide file tree
Showing 7 changed files with 38 additions and 41 deletions.
12 changes: 6 additions & 6 deletions base/linalg/bidiag.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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)
Expand All @@ -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
Expand Down
2 changes: 1 addition & 1 deletion base/linalg/factorization.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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)) :
Expand Down
4 changes: 0 additions & 4 deletions base/linalg/generic.jl
Original file line number Diff line number Diff line change
@@ -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)

Expand Down
2 changes: 1 addition & 1 deletion base/linalg/sparse.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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())
Expand Down
14 changes: 7 additions & 7 deletions base/linalg/special.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand All @@ -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))
Expand Down Expand Up @@ -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)

43 changes: 22 additions & 21 deletions base/linalg/triangular.jl
Original file line number Diff line number Diff line change
@@ -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)
Expand Down Expand Up @@ -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

Expand Down Expand Up @@ -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]))
Expand Down Expand Up @@ -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)
Expand All @@ -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'),
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -733,20 +734,20 @@ 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

## The general promotion methods
### 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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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.")
Expand All @@ -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

2 changes: 1 addition & 1 deletion base/linalg/uniformscaling.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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.λ

Expand Down

0 comments on commit e40c11b

Please sign in to comment.