Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

making Bi/Tri/Sym times Diag output a structured matrix (not sparse) #31889

Merged
merged 12 commits into from
Jul 6, 2019
Merged

making Bi/Tri/Sym times Diag output a structured matrix (not sparse) #31889

Show file tree
Hide file tree
Changes from 11 commits
Commits
Show all changes
12 commits
Select commit Hold shift + click to select a range
15e56c9
making Bi/Tri/Sym times Diag output a structured matrix (not sparse)
mcognetta Apr 30, 2019
aa30721
fix test errors/typos
mcognetta May 1, 2019
440ad8a
Merge branch 'master' of https://github.com/JuliaLang/julia into tri_…
mcognetta May 2, 2019
9250a0d
adding diag*bi/tri/sym and removing unused union
mcognetta May 2, 2019
f4896d8
add mul!(C,::BiTriSym, ::*Diagonal)
dkarrasch May 3, 2019
b59d135
Merge pull request #1 from dkarrasch/patch-1
mcognetta May 3, 2019
801c09d
Merge branch 'master' into tri_times_diag
mcognetta May 18, 2019
ef8019e
Merge remote-tracking branch 'origin' into tri_times_diag
mcognetta Jun 14, 2019
cab2dec
fixing ambiguious methods and adjoint/transpose
mcognetta Jun 14, 2019
91a21bc
space
mcognetta Jun 15, 2019
9c917c6
test error
mcognetta Jun 15, 2019
f4b7f9b
remove transpose
mcognetta Jun 19, 2019
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
87 changes: 78 additions & 9 deletions stdlib/LinearAlgebra/src/bidiag.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -336,7 +336,6 @@ end
const BiTriSym = Union{Bidiagonal,Tridiagonal,SymTridiagonal}
const BiTri = Union{Bidiagonal,Tridiagonal}
mul!(C::AbstractMatrix, A::SymTridiagonal, B::BiTriSym) = A_mul_B_td!(C, A, B)
mul!(C::AbstractMatrix, A::BiTri, B::BiTriSym) = A_mul_B_td!(C, A, B)
mul!(C::AbstractMatrix, A::BiTriSym, B::BiTriSym) = A_mul_B_td!(C, A, B)
mul!(C::AbstractMatrix, A::AbstractTriangular, B::BiTriSym) = A_mul_B_td!(C, A, B)
mul!(C::AbstractMatrix, A::AbstractMatrix, B::BiTriSym) = A_mul_B_td!(C, A, B)
Expand All @@ -347,12 +346,12 @@ mul!(C::AbstractMatrix, A::Adjoint{<:Any,<:AbstractTriangular}, B::BiTriSym) = A
mul!(C::AbstractMatrix, A::Transpose{<:Any,<:AbstractTriangular}, B::BiTriSym) = A_mul_B_td!(C, A, B)
mul!(C::AbstractMatrix, A::Adjoint{<:Any,<:AbstractVecOrMat}, B::BiTriSym) = A_mul_B_td!(C, A, B)
mul!(C::AbstractMatrix, A::Transpose{<:Any,<:AbstractVecOrMat}, B::BiTriSym) = A_mul_B_td!(C, A, B)
mul!(C::AbstractVector, A::BiTri, B::AbstractVector) = A_mul_B_td!(C, A, B)
mul!(C::AbstractMatrix, A::BiTri, B::AbstractVecOrMat) = A_mul_B_td!(C, A, B)
mul!(C::AbstractVecOrMat, A::BiTri, B::AbstractVecOrMat) = A_mul_B_td!(C, A, B)
mul!(C::AbstractMatrix, A::BiTri, B::Transpose{<:Any,<:AbstractVecOrMat}) = A_mul_B_td!(C, A, B) # around bidiag line 330
mul!(C::AbstractMatrix, A::BiTri, B::Adjoint{<:Any,<:AbstractVecOrMat}) = A_mul_B_td!(C, A, B)
mul!(C::AbstractVector, A::BiTri, B::Transpose{<:Any,<:AbstractVecOrMat}) = throw(MethodError(mul!, (C, A, B)))
mul!(C::AbstractVector, A::BiTriSym, B::AbstractVector) = A_mul_B_td!(C, A, B)
mul!(C::AbstractMatrix, A::BiTriSym, B::AbstractVecOrMat) = A_mul_B_td!(C, A, B)
mul!(C::AbstractVecOrMat, A::BiTriSym, B::AbstractVecOrMat) = A_mul_B_td!(C, A, B)
mul!(C::AbstractMatrix, A::BiTriSym, B::Transpose{<:Any,<:AbstractVecOrMat}) = A_mul_B_td!(C, A, B) # around bidiag line 330
mul!(C::AbstractMatrix, A::BiTriSym, B::Adjoint{<:Any,<:AbstractVecOrMat}) = A_mul_B_td!(C, A, B)
mul!(C::AbstractVector, A::BiTriSym, B::Transpose{<:Any,<:AbstractVecOrMat}) = throw(MethodError(mul!, (C, A, B)))

function check_A_mul_B!_sizes(C, A, B)
require_one_based_indexing(C)
Expand Down Expand Up @@ -437,6 +436,39 @@ function A_mul_B_td!(C::AbstractMatrix, A::BiTriSym, B::BiTriSym)
C
end

function A_mul_B_td!(C::AbstractMatrix, A::BiTriSym, B::Diagonal)
Copy link
Member

Choose a reason for hiding this comment

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

Similarly to the comment below: do we need methods for Transposes and Adjoints of B? In line 634, you call A_mul_B_td! with a B::Transpose{<:Any, <:Diagonal}.

Hm, actually, it seems that the exact same code works for B::Union{Diagonal, Adjoint{<:Any, <:Diagonal}, Transpose{<:Any, <:Diagonal}. The elementwise transposition/adjoint is done at getindex. You could define a constant AdjOrTransDiagonal somewhere above?

check_A_mul_B!_sizes(C, A, B)
n = size(A,1)
n <= 3 && return mul!(C, Array(A), Array(B))
fill!(C, zero(eltype(C)))
Al = _diag(A, -1)
Ad = _diag(A, 0)
Au = _diag(A, 1)
Bd = B.diag
@inbounds begin
# first row of C
C[1,1] = A[1,1]*B[1,1]
C[1,2] = A[1,2]*B[2,2]
# second row of C
C[2,1] = A[2,1]*B[1,1]
C[2,2] = A[2,2]*B[2,2]
C[2,3] = A[2,3]*B[3,3]
for j in 3:n-2
C[j, j-1] = Al[j-1]*Bd[j-1]
C[j, j ] = Ad[j ]*Bd[j ]
C[j, j+1] = Au[j ]*Bd[j+1]
end
# row before last of C
C[n-1,n-2] = A[n-1,n-2]*B[n-2,n-2]
C[n-1,n-1] = A[n-1,n-1]*B[n-1,n-1]
C[n-1,n ] = A[n-1, n]*B[n ,n ]
# last row of C
C[n,n-1] = A[n,n-1]*B[n-1,n-1]
C[n,n ] = A[n,n ]*B[n, n ]
end # inbounds
C
end

function A_mul_B_td!(C::AbstractVecOrMat, A::BiTriSym, B::AbstractVecOrMat)
require_one_based_indexing(C)
require_one_based_indexing(B)
Expand Down Expand Up @@ -497,7 +529,39 @@ function A_mul_B_td!(C::AbstractMatrix, A::AbstractMatrix, B::BiTriSym)
C
end

const SpecialMatrix = Union{Bidiagonal,SymTridiagonal,Tridiagonal}
function A_mul_B_td!(C::AbstractMatrix, A::Diagonal, B::BiTriSym)
check_A_mul_B!_sizes(C, A, B)
n = size(A,1)
n <= 3 && return mul!(C, Array(A), Array(B))
fill!(C, zero(eltype(C)))
Ad = A.diag
Bl = _diag(B, -1)
Bd = _diag(B, 0)
Bu = _diag(B, 1)
@inbounds begin
# first row of C
C[1,1] = A[1,1]*B[1,1]
C[1,2] = A[1,1]*B[1,2]
# second row of C
C[2,1] = A[2,2]*B[2,1]
C[2,2] = A[2,2]*B[2,2]
C[2,3] = A[2,2]*B[2,3]
for j in 3:n-2
Ajj = Ad[j]
C[j, j-1] = Ajj*Bl[j-1]
C[j, j ] = Ajj*Bd[j]
C[j, j+1] = Ajj*Bu[j]
end
# row before last of C
C[n-1,n-2] = A[n-1,n-1]*B[n-1,n-2]
C[n-1,n-1] = A[n-1,n-1]*B[n-1,n-1]
C[n-1,n ] = A[n-1,n-1]*B[n-1,n ]
# last row of C
C[n,n-1] = A[n,n]*B[n,n-1]
C[n,n ] = A[n,n]*B[n,n ]
end # inbounds
C
end

function *(A::AbstractTriangular, B::Union{SymTridiagonal, Tridiagonal})
TS = promote_op(matprod, eltype(A), eltype(B))
Expand Down Expand Up @@ -545,7 +609,7 @@ function *(A::Bidiagonal, B::LowerTriangular)
end
end

function *(A::Bidiagonal, B::Diagonal)
function *(A::BiTri, B::Diagonal)
TS = promote_op(matprod, eltype(A), eltype(B))
A_mul_B_td!(similar(A, TS), A, B)
end
Expand All @@ -565,6 +629,11 @@ function *(A::SymTridiagonal, B::Diagonal)
A_mul_B_td!(Tridiagonal(zeros(TS, size(A, 1)-1), zeros(TS, size(A, 1)), zeros(TS, size(A, 1)-1)), A, B)
end

function *(A::BiTri, B::Transpose{<:Any, <:Diagonal})
TS = promote_op(matprod, eltype(A), eltype(B))
A_mul_B_td!(similar(A, TS), A, B)
end

Copy link
Member

Choose a reason for hiding this comment

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

What about B::Adjoint{<:Any, <:Diagonal}? And what about A::SymTridiagonal? Are these handled by other * methods?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

The thing is that even if we add all these Transpose/Adjoint methods, they will be slow due to similar(Transpose(D), Float64) returning a SparseMatrix. This is what I was alluding to in #31889 (comment).

Personally, I think that similar for structured matrices should always return the same type (unless the dims are not square), but this seems to be a design decision that a lot of thought went into.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Perhaps this should be removed and put into another PR where we handle the Adjoint/Transpose case?

#Generic multiplication
*(A::Bidiagonal{T}, B::AbstractVector{T}) where {T} = *(Array(A), B)
*(adjA::Adjoint{<:Any,<:Bidiagonal{T}}, B::AbstractVector{T}) where {T} = *(adjoint(Array(adjA.parent)), B)
Expand Down
32 changes: 32 additions & 0 deletions stdlib/LinearAlgebra/test/bidiag.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -304,6 +304,38 @@ Random.seed!(1)
C = Matrix{elty}(undef, n, n)
Dia = Diagonal(T.dv)
@test mul!(C, Dia, T) ≈ Array(Dia)*Array(T)

# Issue #31870
# Bi/Tri/Sym times Diagonal
Diag = Diagonal(rand(elty, 10))
BidiagU = Bidiagonal(rand(elty, 10), rand(elty, 9), 'U')
BidiagL = Bidiagonal(rand(elty, 10), rand(elty, 9), 'L')
Tridiag = Tridiagonal(rand(elty, 9), rand(elty, 10), rand(elty, 9))
SymTri = SymTridiagonal(rand(elty, 10), rand(elty, 9))

mats = [Diag, BidiagU, BidiagL, Tridiag, SymTri]
for a in mats
for b in mats
for lwrap in (Transpose, Adjoint)
for rwrap in (Transpose, Adjoint)
@test a*b ≈ Matrix(a)*Matrix(b)
@test lwrap(a)*b ≈ Matrix(lwrap(a))*Matrix(b)
@test a*rwrap(b) ≈ Matrix(a)*Matrix(rwrap(b))
@test lwrap(a)*rwrap(b) ≈ Matrix(lwrap(a))*Matrix(rwrap(b))
end
end
end
end

@test typeof(BidiagU*Diag) <: Bidiagonal
@test typeof(BidiagL*Diag) <: Bidiagonal
@test typeof(Tridiag*Diag) <: Tridiagonal
@test typeof(SymTri*Diag) <: Tridiagonal

@test typeof(BidiagU*Diag) <: Bidiagonal
@test typeof(Diag*BidiagL) <: Bidiagonal
@test typeof(Diag*Tridiag) <: Tridiagonal
@test typeof(Diag*SymTri) <: Tridiagonal
end

@test inv(T)*Tfull ≈ Matrix(I, n, n)
Expand Down