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.

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

lazier, less-jazzy linalg internals #24969

Merged
merged 65 commits into from
Dec 12, 2017
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
65 commits
Select commit Hold shift + click to select a range
b0429ff
Basic definitions and tests for Adjoint and Transpose types.
Sacha0 Dec 4, 2017
06bf431
Replace A[ct]_(mul|ldiv|rdiv)_B[ct][!] defs in base/linalg/bidiag wit…
Sacha0 Dec 4, 2017
2208762
Replace A[ct]_(mul|ldiv|rdiv)_B[ct][!] defs in base/linalg/tridiag.jl…
Sacha0 Dec 4, 2017
39cc5cc
Replace A[ct]_(mul|ldiv|rdiv)_B[ct][!] defs in base/linalg/diagonal.j…
Sacha0 Dec 4, 2017
badadbf
Replace A[ct]_(mul|ldiv|rdiv)_B[ct][!] defs in base/linalg/special.jl…
Sacha0 Dec 4, 2017
03057a9
Replace A[ct]_(mul|ldiv|rdiv)_B[ct][!] defs in base/linalg/bunchkaufm…
Sacha0 Dec 4, 2017
89cc4e7
Replace A[ct]_(mul|ldiv|rdiv)_B[ct][!] defs in base/linalg/cholesky.j…
Sacha0 Dec 4, 2017
c138fdc
Replace A[ct]_(mul|ldiv|rdiv)_B[ct][!] defs in base/linalg/factorizat…
Sacha0 Dec 5, 2017
4642c3e
Replace A[ct]_(mul|ldiv|rdiv)_B[ct][!] defs in base/linalg/hessenberg…
Sacha0 Dec 5, 2017
e7f2125
Replace A[ct]_(mul|ldiv|rdiv)_B[ct][!] defs in base/linalg/ldlt.jl wi…
Sacha0 Dec 5, 2017
5c45fef
Replace A[ct]_(mul|ldiv|rdiv)_B[ct][!] defs in base/linalg/svd.jl wit…
Sacha0 Dec 5, 2017
f5f77a5
Replace A[ct]_(mul|ldiv|rdiv)_B[ct][!] defs in base/linalg/symmetric.…
Sacha0 Dec 5, 2017
6fdb961
Replace A[ct]_(mul|ldiv|rdiv)_B[ct][!] defs in base/linalg/lu.jl with…
Sacha0 Dec 5, 2017
79a2ce2
Replace A[ct]_(mul|ldiv|rdiv)_B[ct][!] defs in base/linalg/lq.jl with…
Sacha0 Dec 5, 2017
096fc67
Replace A[ct]_(mul|ldiv|rdiv)_B[ct][!] defs in base/linalg/qr.jl with…
Sacha0 Dec 5, 2017
2a29a50
Replace A[ct]_(mul|ldiv|rdiv)_B[ct][!] defs in base/linalg/matmul.jl …
Sacha0 Dec 5, 2017
e41efab
Replace A[ct]_(mul|ldiv|rdiv)_B[ct][!] defs in base/linalg/triangular…
Sacha0 Dec 6, 2017
57183d2
Replace A[ct]_(mul|ldiv|rdiv)_B[ct][!] defs in base/sparse/linalg.jl …
Sacha0 Dec 6, 2017
ae76f0d
Replace A[ct]_(mul|ldiv|rdiv)_B[ct][!] defs in base/sparse/sparsevect…
Sacha0 Dec 6, 2017
2edbb1c
Replace A[ct]_(mul|ldiv|rdiv)_B[ct][!] defs in base/linalg/rowvector.…
Sacha0 Dec 6, 2017
5163363
Replace A[ct]_(mul|ldiv|rdiv)_B[ct][!] defs in base/linalg/givens.jl …
Sacha0 Dec 7, 2017
e524a32
Replace A[ct]_(mul|ldiv|rdiv)_B[ct][!] defs in stdlib/IterativeEigenS…
Sacha0 Dec 9, 2017
8d20c31
Replace A[ct]_(mul|ldiv|rdiv)_B[ct][!] defs in stdlib/SuiteSparse/CHO…
Sacha0 Dec 10, 2017
0d4f196
Replace A[ct]_(mul|ldiv|rdiv)_B[ct][!] defs in stdlib/SuiteSparse/UMF…
Sacha0 Dec 10, 2017
d062298
Replace A[ct]_(mul|ldiv|rdiv)_B[ct][!] defs in stdlib/SuiteSparse/SPQ…
Sacha0 Dec 10, 2017
ce5ee65
Disambiguate *, /, and \ methods.
Sacha0 Dec 9, 2017
6fb85f0
Disambiguate mul!, ldiv!, and rdiv! methods.
Sacha0 Dec 10, 2017
e80697d
Replace A[ct]_(mul|ldiv|rdiv)_B[ct][!] defs in base/operators.jl with…
Sacha0 Dec 8, 2017
3f2d810
Rewrite A[ct]_(mul|ldiv|rdiv)_B[ct][!] calls in base/deprecated.jl as…
Sacha0 Dec 11, 2017
233175d
Rewrite A[ct]_(mul|ldiv|rdiv)_B[ct][!] calls in base/linalg/bidiag.jl…
Sacha0 Dec 11, 2017
c1aa72c
Rewrite A[ct]_(mul|ldiv|rdiv)_B[ct][!] calls in base/linalg/bunchkauf…
Sacha0 Dec 11, 2017
ecd0253
Rewrite A[ct]_(mul|ldiv|rdiv)_B[ct][!] calls in base/linalg/cholesky.…
Sacha0 Dec 11, 2017
fc0accd
Rewrite A[ct]_(mul|ldiv|rdiv)_B[ct][!] calls in base/linalg/dense.jl …
Sacha0 Dec 11, 2017
5673e1a
Rewrite A[ct]_(mul|ldiv|rdiv)_B[ct][!] calls in base/linalg/diagonal.…
Sacha0 Dec 11, 2017
c98439d
Rewrite A[ct]_(mul|ldiv|rdiv)_B[ct][!] calls in base/linalg/factoriza…
Sacha0 Dec 11, 2017
a12f57e
Rewrite A[ct]_(mul|ldiv|rdiv)_B[ct][!] calls in base/linalg/generic.j…
Sacha0 Dec 11, 2017
16c5ce7
Rewrite A[ct]_(mul|ldiv|rdiv)_B[ct][!] calls in base/linalg/givens.jl…
Sacha0 Dec 11, 2017
9e9afcf
Rewrite A[ct]_(mul|ldiv|rdiv)_B[ct][!] calls in base/linalg/hessenber…
Sacha0 Dec 11, 2017
aa56327
Rewrite A[ct]_(mul|ldiv|rdiv)_B[ct][!] calls in base/linalg/lq.jl as …
Sacha0 Dec 11, 2017
2883770
Rewrite A[ct]_(mul|ldiv|rdiv)_B[ct][!] calls in base/linalg/lu.jl as …
Sacha0 Dec 11, 2017
5226436
Rewrite A[ct]_(mul|ldiv|rdiv)_B[ct][!] calls in base/linalg/matmul.jl…
Sacha0 Dec 11, 2017
e1db456
Rewrite A[ct]_(mul|ldiv|rdiv)_B[ct][!] calls in base/linalg/qr.jl as …
Sacha0 Dec 11, 2017
d2a4b16
Rewrite A[ct]_(mul|ldiv|rdiv)_B[ct][!] calls in base/linalg/special.j…
Sacha0 Dec 11, 2017
591d226
Rewrite A[ct]_(mul|ldiv|rdiv)_B[ct][!] calls in base/sparse/linalg.jl…
Sacha0 Dec 11, 2017
be73121
Rewrite A[ct]_(mul|ldiv|rdiv)_B[ct][!] calls in base/statistics.jl as…
Sacha0 Dec 11, 2017
bdba40b
Rewrite A[ct]_(mul|ldiv|rdiv)_B[ct][!] calls in base/sparse/sparsevec…
Sacha0 Dec 11, 2017
9b387f9
Rewrite A[ct]_(mul|ldiv|rdiv)_B[ct][!] calls in base/linalg/triangula…
Sacha0 Dec 11, 2017
5147f52
Rewrite A[ct]_(mul|ldiv|rdiv)_B[ct][!] calls in base/linalg/rowvector…
Sacha0 Dec 11, 2017
364a144
Rewrite A[ct]_(mul|ldiv|rdiv)_B[ct][!] calls in test/linalg/diagonal.…
Sacha0 Dec 11, 2017
8626205
Rewrite A[ct]_(mul|ldiv|rdiv)_B[ct][!] calls in test/linalg/bidiag.jl…
Sacha0 Dec 11, 2017
116b951
Rewrite A[ct]_(mul|ldiv|rdiv)_B[ct][!] calls in test/linalg/givens.jl…
Sacha0 Dec 11, 2017
8b773fe
Rewrite A[ct]_(mul|ldiv|rdiv)_B[ct][!] calls in test/linalg/lq.jl as …
Sacha0 Dec 11, 2017
613f320
Rewrite A[ct]_(mul|ldiv|rdiv)_B[ct][!] calls in test/linalg/lu.jl as …
Sacha0 Dec 11, 2017
306d1dd
Rewrite A[ct]_(mul|ldiv|rdiv)_B[ct][!] calls in test/linalg/matmul.jl…
Sacha0 Dec 11, 2017
cf4988e
Rewrite A[ct]_(mul|ldiv|rdiv)_B[ct][!] calls in test/linalg/qr.jl as …
Sacha0 Dec 11, 2017
5e9a327
Rewrite A[ct]_(mul|ldiv|rdiv)_B[ct][!] calls in test/linalg/special.j…
Sacha0 Dec 11, 2017
6e0f0ae
Rewrite A[ct]_(mul|ldiv|rdiv)_B[ct][!] calls in test/linalg/symmetric…
Sacha0 Dec 11, 2017
df92f58
Rewrite A[ct]_(mul|ldiv|rdiv)_B[ct][!] calls in test/linalg/triangula…
Sacha0 Dec 11, 2017
25f36af
Rewrite A[ct]_(mul|ldiv|rdiv)_B[ct][!] calls in test/linalg/tridiag.j…
Sacha0 Dec 11, 2017
061c6bb
Rewrite A[ct]_(mul|ldiv|rdiv)_B[ct][!] calls in test/linalg/uniformsc…
Sacha0 Dec 11, 2017
b306292
Rewrite A[ct]_(mul|ldiv|rdiv)_B[ct][!] calls in test/perf/* as *, /, …
Sacha0 Dec 11, 2017
2b91b2f
Rewrite A[ct]_(mul|ldiv|rdiv)_B[ct][!] calls in test/sparse/sparse.jl…
Sacha0 Dec 11, 2017
fa97d9c
Rewrite A[ct]_(mul|ldiv|rdiv)_B[ct][!] calls in test/sparse/sparsevec…
Sacha0 Dec 11, 2017
a11cbe5
Rewrite A[ct]_(mul|ldiv|rdiv)_B[ct][!] calls in stdlib/IterativeEigen…
Sacha0 Dec 12, 2017
a05e85f
Rewrite A[ct]_(mul|ldiv|rdiv)_B[ct][!] calls in stdlib/SuiteSparse as…
Sacha0 Dec 12, 2017
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
794 changes: 790 additions & 4 deletions base/deprecated.jl

Large diffs are not rendered by default.

120 changes: 120 additions & 0 deletions base/linalg/adjtrans.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,120 @@
# This file is a part of Julia. License is MIT: https://julialang.org/license

using Base: @pure, @propagate_inbounds, _return_type, _default_type, _isleaftype, @_inline_meta
import Base: length, size, indices, IndexStyle, getindex, setindex!, parent, vec, convert, similar

### basic definitions (types, aliases, constructors, abstractarray interface, sundry similar)

# note that Adjoint and Transpose must be able to wrap not only vectors and matrices
# but also factorizations, rotations, and other linear algebra objects, including
# user-defined such objects. so do not restrict the wrapped type.
struct Adjoint{T,S} <: AbstractMatrix{T}
parent::S
function Adjoint{T,S}(A::S) where {T,S}
checkeltype(Adjoint, T, eltype(A))
new(A)
end
end
struct Transpose{T,S} <: AbstractMatrix{T}
parent::S
function Transpose{T,S}(A::S) where {T,S}
checkeltype(Transpose, T, eltype(A))
new(A)
end
end

@pure function checkeltype(::Type{Transform}, ::Type{ResultEltype}, ::Type{ParentEltype}) where {Transform, ResultEltype, ParentEltype}
Copy link
Member

Choose a reason for hiding this comment

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

Why does this need to be @pure

Copy link
Member Author

Choose a reason for hiding this comment

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

This definition was preexisting for RowVectors :).

Copy link
Member

Choose a reason for hiding this comment

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

This is very strongly not a pure function and terrible things will happen.

Copy link
Member Author

Choose a reason for hiding this comment

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

As above, this definition was preexisting for RowVectors :).

if ResultEltype !== transformtype(Transform, ParentEltype)
error(string("Element type mismatch. Tried to create an `$Transform{$ResultEltype}` ",
"from an object with eltype `$ParentEltype`, but the element type of the ",
"`$Transform` of an object with eltype `$ParentEltype` must be ",
"`$(transformtype(Transform, ParentEltype))`"))
end
return nothing
end
function transformtype(::Type{O}, ::Type{S}) where {O,S}
# similar to promote_op(::Any, ::Type)
@_inline_meta
T = _return_type(O, Tuple{_default_type(S)})
_isleaftype(S) && return _isleaftype(T) ? T : Any
return typejoin(S, T)
end

# basic outer constructors
Adjoint(A) = Adjoint{transformtype(Adjoint,eltype(A)),typeof(A)}(A)
Transpose(A) = Transpose{transformtype(Transpose,eltype(A)),typeof(A)}(A)

# numbers are the end of the line
Adjoint(x::Number) = adjoint(x)
Transpose(x::Number) = transpose(x)

# unwrapping constructors
# perhaps slightly odd, but necessary (at least till adjoint and transpose names are free)
Adjoint(A::Adjoint) = A.parent
Transpose(A::Transpose) = A.parent

# some aliases for internal convenience use
const AdjOrTrans{T,S} = Union{Adjoint{T,S},Transpose{T,S}} where {T,S}
const AdjOrTransAbsVec{T} = AdjOrTrans{T,<:AbstractVector}
const AdjOrTransAbsMat{T} = AdjOrTrans{T,<:AbstractMatrix}

# for internal use below
wrappertype(A::Adjoint) = Adjoint
wrappertype(A::Transpose) = Transpose
wrappertype(::Type{<:Adjoint}) = Adjoint
wrappertype(::Type{<:Transpose}) = Transpose

# AbstractArray interface, basic definitions
length(A::AdjOrTrans) = length(A.parent)
size(v::AdjOrTransAbsVec) = (1, length(v.parent))
size(A::AdjOrTransAbsMat) = reverse(size(A.parent))
indices(v::AdjOrTransAbsVec) = (Base.OneTo(1), indices(v.parent)...)
indices(A::AdjOrTransAbsMat) = reverse(indices(A.parent))
IndexStyle(::Type{<:AdjOrTransAbsVec}) = IndexLinear()
IndexStyle(::Type{<:AdjOrTransAbsMat}) = IndexCartesian()
@propagate_inbounds getindex(v::AdjOrTransAbsVec, i::Int) = wrappertype(v)(v.parent[i])
@propagate_inbounds getindex(A::AdjOrTransAbsMat, i::Int, j::Int) = wrappertype(A)(A.parent[j, i])
@propagate_inbounds setindex!(v::AdjOrTransAbsVec, x, i::Int) = (setindex!(v.parent, wrappertype(v)(x), i); v)
@propagate_inbounds setindex!(A::AdjOrTransAbsMat, x, i::Int, j::Int) = (setindex!(A.parent, wrappertype(A)(x), j, i); A)
# AbstractArray interface, additional definitions to retain wrapper over vectors where appropriate
@propagate_inbounds getindex(v::AdjOrTransAbsVec, ::Colon, is::AbstractArray{Int}) = wrappertype(v)(v.parent[is])
@propagate_inbounds getindex(v::AdjOrTransAbsVec, ::Colon, ::Colon) = wrappertype(v)(v.parent[:])

# conversion of underlying storage
convert(::Type{Adjoint{T,S}}, A::Adjoint) where {T,S} = Adjoint{T,S}(convert(S, A.parent))
convert(::Type{Transpose{T,S}}, A::Transpose) where {T,S} = Transpose{T,S}(convert(S, A.parent))

# for vectors, the semantics of the wrapped and unwrapped types differ
# so attempt to maintain both the parent and wrapper type insofar as possible
similar(A::AdjOrTransAbsVec) = wrappertype(A)(similar(A.parent))
similar(A::AdjOrTransAbsVec, ::Type{T}) where {T} = wrappertype(A)(similar(A.parent, transformtype(wrappertype(A), T)))
# for matrices, the semantics of the wrapped and unwrapped types are generally the same
# and as you are allocating with similar anyway, you might as well get something unwrapped
similar(A::AdjOrTrans) = similar(A.parent, eltype(A), size(A))
similar(A::AdjOrTrans, ::Type{T}) where {T} = similar(A.parent, T, size(A))
similar(A::AdjOrTrans, ::Type{T}, dims::Dims{N}) where {T,N} = similar(A.parent, T, dims)

# sundry basic definitions
parent(A::AdjOrTrans) = A.parent
vec(v::AdjOrTransAbsVec) = v.parent


### linear algebra

# definitions necessary for test/linalg/rowvector.jl to pass
# should be cleaned up / revised as necessary in the future
/(A::Transpose{<:Any,<:Vector}, B::Matrix) = /(transpose(A.parent), B)
/(A::Transpose{<:Any,<:Vector}, B::Transpose{<:Any,<:Matrix}) = /(transpose(A.parent), B)
*(A::Adjoint{<:Any,<:Matrix}, B::Adjoint{<:Any,<:Vector}) = *(adjoint(A.parent), adjoint(B.parent))


# dismabiguation methods
*(A::Transpose{<:Any,<:AbstractVector}, B::Adjoint{<:Any,<:AbstractVector}) = transpose(A.parent) * B
*(A::Transpose{<:Any,<:AbstractVector}, B::Adjoint{<:Any,<:AbstractMatrix}) = transpose(A.parent) * B
*(A::Transpose{<:Any,<:AbstractMatrix}, B::Adjoint{<:Any,<:AbstractVector}) = A * adjoint(B.parent)
*(A::Transpose{<:Any,<:AbstractMatrix}, B::Adjoint{<:Any,<:AbstractMatrix}) = transpose(A.parent) * B
*(A::Adjoint{<:Any,<:AbstractVector}, B::Transpose{<:Any,<:AbstractVector}) = adjoint(A.parent) * B
*(A::Adjoint{<:Any,<:AbstractVector}, B::Transpose{<:Any,<:AbstractMatrix}) = adjoint(A.parent) * B
*(A::Adjoint{<:Any,<:AbstractMatrix}, B::Adjoint{<:Any,<:AbstractVector}) = A * adjoint(B.parent)
*(A::Adjoint{<:Any,<:AbstractMatrix}, B::Transpose{<:Any,<:AbstractVector}) = A * transpose(B.parent)
*(A::Adjoint{<:Any,<:AbstractMatrix}, B::Transpose{<:Any,<:AbstractMatrix}) = adjoint(A.parent) * B
139 changes: 87 additions & 52 deletions base/linalg/bidiag.jl
Original file line number Diff line number Diff line change
Expand Up @@ -325,25 +325,33 @@ end

const BiTriSym = Union{Bidiagonal,Tridiagonal,SymTridiagonal}
const BiTri = Union{Bidiagonal,Tridiagonal}
A_mul_B!(C::AbstractMatrix, A::SymTridiagonal, B::BiTriSym) = A_mul_B_td!(C, A, B)
A_mul_B!(C::AbstractMatrix, A::BiTri, B::BiTriSym) = A_mul_B_td!(C, A, B)
A_mul_B!(C::AbstractMatrix, A::BiTriSym, B::BiTriSym) = A_mul_B_td!(C, A, B)
A_mul_B!(C::AbstractMatrix, A::AbstractTriangular, B::BiTriSym) = A_mul_B_td!(C, A, B)
A_mul_B!(C::AbstractMatrix, A::AbstractMatrix, B::BiTriSym) = A_mul_B_td!(C, A, B)
A_mul_B!(C::AbstractMatrix, A::Diagonal, B::BiTriSym) = A_mul_B_td!(C, A, B)
A_mul_B!(C::AbstractVector, A::BiTri, B::AbstractVector) = A_mul_B_td!(C, A, B)
A_mul_B!(C::AbstractMatrix, A::BiTri, B::AbstractVecOrMat) = A_mul_B_td!(C, A, B)
A_mul_B!(C::AbstractVecOrMat, A::BiTri, B::AbstractVecOrMat) = A_mul_B_td!(C, A, B)

\(::Diagonal, ::RowVector) = throw(DimensionMismatch("Cannot left-divide matrix by transposed vector"))
\(::Bidiagonal, ::RowVector) = throw(DimensionMismatch("Cannot left-divide matrix by transposed vector"))
\(::Bidiagonal{<:Number}, ::RowVector{<:Number}) = throw(DimensionMismatch("Cannot left-divide matrix by transposed vector"))

At_ldiv_B(::Bidiagonal, ::RowVector) = throw(DimensionMismatch("Cannot left-divide matrix by transposed vector"))
At_ldiv_B(::Bidiagonal{<:Number}, ::RowVector{<:Number}) = throw(DimensionMismatch("Cannot left-divide matrix by transposed vector"))

Ac_ldiv_B(::Bidiagonal, ::RowVector) = throw(DimensionMismatch("Cannot left-divide matrix by transposed vector"))
Ac_ldiv_B(::Bidiagonal{<:Number}, ::RowVector{<:Number}) = throw(DimensionMismatch("Cannot left-divide matrix by transposed vector"))
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)
mul!(C::AbstractMatrix, A::Diagonal, B::BiTriSym) = A_mul_B_td!(C, A, B)
mul!(C::AbstractMatrix, A::Adjoint{<:Any,<:Diagonal}, B::BiTriSym) = A_mul_B_td!(C, A, B)
mul!(C::AbstractMatrix, A::Transpose{<:Any,<:Diagonal}, B::BiTriSym) = A_mul_B_td!(C, A, B)
mul!(C::AbstractMatrix, A::Adjoint{<:Any,<:AbstractTriangular}, B::BiTriSym) = A_mul_B_td!(C, A, B)
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)))

\(::Diagonal, ::RowVector) = _mat_ldiv_rowvec_error()
\(::Bidiagonal, ::RowVector) = _mat_ldiv_rowvec_error()
\(::Bidiagonal{<:Number}, ::RowVector{<:Number}) = _mat_ldiv_rowvec_error()
\(::Adjoint{<:Any,<:Bidiagonal}, ::RowVector) = _mat_ldiv_rowvec_error()
\(::Transpose{<:Any,<:Bidiagonal}, ::RowVector) = _mat_ldiv_rowvec_error()
\(::Adjoint{<:Number,<:Bidiagonal{<:Number}}, ::RowVector{<:Number}) = _mat_ldiv_rowvec_error()
\(::Transpose{<:Number,<:Bidiagonal{<:Number}}, ::RowVector{<:Number}) = _mat_ldiv_rowvec_error()
_mat_ldiv_rowvec_error() = throw(DimensionMismatch("Cannot left-divide matrix by transposed vector"))

function check_A_mul_B!_sizes(C, A, B)
nA, mA = size(A)
Expand Down Expand Up @@ -375,7 +383,7 @@ end
function A_mul_B_td!(C::AbstractMatrix, A::BiTriSym, B::BiTriSym)
check_A_mul_B!_sizes(C, A, B)
n = size(A,1)
n <= 3 && return A_mul_B!(C, Array(A), Array(B))
n <= 3 && return mul!(C, Array(A), Array(B))
fill!(C, zero(eltype(C)))
Al = _diag(A, -1)
Ad = _diag(A, 0)
Expand Down Expand Up @@ -434,7 +442,7 @@ function A_mul_B_td!(C::AbstractVecOrMat, A::BiTriSym, B::AbstractVecOrMat)
if size(C,2) != nB
throw(DimensionMismatch("A has second dimension $nA, B has $(size(B,2)), C has $(size(C,2)) but all must match"))
end
nA <= 3 && return A_mul_B!(C, Array(A), Array(B))
nA <= 3 && return mul!(C, Array(A), Array(B))
l = _diag(A, -1)
d = _diag(A, 0)
u = _diag(A, 1)
Expand All @@ -455,7 +463,7 @@ end
function A_mul_B_td!(C::AbstractMatrix, A::AbstractMatrix, B::BiTriSym)
check_A_mul_B!_sizes(C, A, B)
n = size(A,1)
n <= 3 && return A_mul_B!(C, Array(A), Array(B))
n <= 3 && return mul!(C, Array(A), Array(B))
m = size(B,2)
Bl = _diag(B, -1)
Bd = _diag(B, 0)
Expand Down Expand Up @@ -489,15 +497,17 @@ const SpecialMatrix = Union{Bidiagonal,SymTridiagonal,Tridiagonal}
*(A::SpecialMatrix, B::SpecialMatrix) = Array(A) * Array(B)

#Generic multiplication
for func in (:*, :Ac_mul_B, :A_mul_Bc, :/, :A_rdiv_Bc)
@eval ($func)(A::Bidiagonal{T}, B::AbstractVector{T}) where {T} = ($func)(Array(A), B)
end
*(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)
*(A::Bidiagonal{T}, adjB::Adjoint{<:Any,<:AbstractVector{T}}) where {T} = *(Array(A), Adjoint(adjB.parent))
/(A::Bidiagonal{T}, B::AbstractVector{T}) where {T} = /(Array(A), B)
/(A::Bidiagonal{T}, adjB::Adjoint{<:Any,<:AbstractVector{T}}) where {T} = /(Array(A), Adjoint(adjB.parent))

#Linear solvers
A_ldiv_B!(A::Union{Bidiagonal, AbstractTriangular}, b::AbstractVector) = naivesub!(A, b)
At_ldiv_B!(A::Bidiagonal, b::AbstractVector) = A_ldiv_B!(transpose(A), b)
Ac_ldiv_B!(A::Bidiagonal, b::AbstractVector) = A_ldiv_B!(adjoint(A), b)
function A_ldiv_B!(A::Union{Bidiagonal,AbstractTriangular}, B::AbstractMatrix)
ldiv!(A::Union{Bidiagonal, AbstractTriangular}, b::AbstractVector) = naivesub!(A, b)
ldiv!(transA::Transpose{<:Any,<:Bidiagonal}, b::AbstractVector) = ldiv!(transpose(transA.parent), b)
ldiv!(adjA::Adjoint{<:Any,<:Bidiagonal}, b::AbstractVector) = ldiv!(adjoint(adjA.parent), b)
function ldiv!(A::Union{Bidiagonal,AbstractTriangular}, B::AbstractMatrix)
nA,mA = size(A)
tmp = similar(B,size(B,1))
n = size(B, 1)
Expand All @@ -506,26 +516,40 @@ function A_ldiv_B!(A::Union{Bidiagonal,AbstractTriangular}, B::AbstractMatrix)
end
for i = 1:size(B,2)
copy!(tmp, 1, B, (i - 1)*n + 1, n)
A_ldiv_B!(A, tmp)
ldiv!(A, tmp)
copy!(B, (i - 1)*n + 1, tmp, 1, n) # Modify this when array view are implemented.
end
B
end
for func in (:Ac_ldiv_B!, :At_ldiv_B!)
@eval function ($func)(A::Union{Bidiagonal,AbstractTriangular}, B::AbstractMatrix)
nA,mA = size(A)
tmp = similar(B,size(B,1))
n = size(B, 1)
if mA != n
throw(DimensionMismatch("size of A' is ($mA,$nA), corresponding dimension of B is $n"))
end
for i = 1:size(B,2)
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
function ldiv!(adjA::Adjoint{<:Any,<:Union{Bidiagonal,AbstractTriangular}}, B::AbstractMatrix)
A = adjA.parent
nA,mA = size(A)
tmp = similar(B,size(B,1))
n = size(B, 1)
if mA != n
throw(DimensionMismatch("size of A' is ($mA,$nA), corresponding dimension of B is $n"))
end
for i = 1:size(B,2)
copy!(tmp, 1, B, (i - 1)*n + 1, n)
ldiv!(Adjoint(A), tmp)
copy!(B, (i - 1)*n + 1, tmp, 1, n) # Modify this when array view are implemented.
end
B
end
function ldiv!(transA::Transpose{<:Any,<:Union{Bidiagonal,AbstractTriangular}}, B::AbstractMatrix)
A = transA.parent
nA,mA = size(A)
tmp = similar(B,size(B,1))
n = size(B, 1)
if mA != n
throw(DimensionMismatch("size of A' is ($mA,$nA), corresponding dimension of B is $n"))
end
for i = 1:size(B,2)
copy!(tmp, 1, B, (i - 1)*n + 1, n)
ldiv!(Transpose(A), tmp)
copy!(B, (i - 1)*n + 1, tmp, 1, n) # Modify this when array view are implemented.
end
B
end
#Generic solver using naive substitution
function naivesub!(A::Bidiagonal{T}, b::AbstractVector, x::AbstractVector = b) where T
Expand All @@ -550,15 +574,26 @@ function naivesub!(A::Bidiagonal{T}, b::AbstractVector, x::AbstractVector = b) w
end

### Generic promotion methods and fallbacks
for (f,g) in ((:\, :A_ldiv_B!), (:At_ldiv_B, :At_ldiv_B!), (:Ac_ldiv_B, :Ac_ldiv_B!))
@eval begin
function ($f)(A::Bidiagonal{TA}, B::AbstractVecOrMat{TB}) where {TA<:Number,TB<:Number}
TAB = typeof((zero(TA)*zero(TB) + zero(TA)*zero(TB))/one(TA))
($g)(convert(AbstractArray{TAB}, A), copy_oftype(B, TAB))
end
($f)(A::Bidiagonal, B::AbstractVecOrMat) = ($g)(A, copy(B))
end
function \(A::Bidiagonal{<:Number}, B::AbstractVecOrMat{<:Number})
TA, TB = eltype(A), eltype(B)
TAB = typeof((zero(TA)*zero(TB) + zero(TA)*zero(TB))/one(TA))
ldiv!(convert(AbstractArray{TAB}, A), copy_oftype(B, TAB))
end
\(A::Bidiagonal, B::AbstractVecOrMat) = ldiv!(A, copy(B))
function \(transA::Transpose{<:Number,<:Bidiagonal{<:Number}}, B::AbstractVecOrMat{<:Number})
A = transA.parent
TA, TB = eltype(A), eltype(B)
TAB = typeof((zero(TA)*zero(TB) + zero(TA)*zero(TB))/one(TA))
ldiv!(Transpose(convert(AbstractArray{TAB}, A)), copy_oftype(B, TAB))
end
\(transA::Transpose{<:Any,<:Bidiagonal}, B::AbstractVecOrMat) = ldiv!(Transpose(transA.parent), copy(B))
function \(adjA::Adjoint{<:Number,<:Bidiagonal{<:Number}}, B::AbstractVecOrMat{<:Number})
A = adjA.parent
TA, TB = eltype(A), eltype(B)
TAB = typeof((zero(TA)*zero(TB) + zero(TA)*zero(TB))/one(TA))
ldiv!(Adjoint(convert(AbstractArray{TAB}, A)), copy_oftype(B, TAB))
end
\(adjA::Adjoint{<:Any,<:Bidiagonal}, B::AbstractVecOrMat) = ldiv!(Adjoint(adjA.parent), copy(B))

factorize(A::Bidiagonal) = A

Expand Down
8 changes: 4 additions & 4 deletions base/linalg/bunchkaufman.jl
Original file line number Diff line number Diff line change
Expand Up @@ -254,7 +254,7 @@ function inv(B::BunchKaufman{<:BlasComplex})
end
end

function A_ldiv_B!(B::BunchKaufman{T}, R::StridedVecOrMat{T}) where T<:BlasReal
function ldiv!(B::BunchKaufman{T}, R::StridedVecOrMat{T}) where T<:BlasReal
if !issuccess(B)
throw(SingularException(B.info))
end
Expand All @@ -265,7 +265,7 @@ function A_ldiv_B!(B::BunchKaufman{T}, R::StridedVecOrMat{T}) where T<:BlasReal
LAPACK.sytrs!(B.uplo, B.LD, B.ipiv, R)
end
end
function A_ldiv_B!(B::BunchKaufman{T}, R::StridedVecOrMat{T}) where T<:BlasComplex
function ldiv!(B::BunchKaufman{T}, R::StridedVecOrMat{T}) where T<:BlasComplex
if !issuccess(B)
throw(SingularException(B.info))
end
Expand All @@ -285,9 +285,9 @@ function A_ldiv_B!(B::BunchKaufman{T}, R::StridedVecOrMat{T}) where T<:BlasCompl
end
end
# There is no fallback solver for Bunch-Kaufman so we'll have to promote to same element type
function A_ldiv_B!(B::BunchKaufman{T}, R::StridedVecOrMat{S}) where {T,S}
function ldiv!(B::BunchKaufman{T}, R::StridedVecOrMat{S}) where {T,S}
TS = promote_type(T,S)
return A_ldiv_B!(convert(BunchKaufman{TS}, B), convert(AbstractArray{TS}, R))
return ldiv!(convert(BunchKaufman{TS}, B), convert(AbstractArray{TS}, R))
end

function logabsdet(F::BunchKaufman)
Expand Down
Loading