Skip to content

Loosen signature in triangular solver from Strided- to AbstractMatrix. #16205

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

Merged
merged 2 commits into from
May 18, 2016
Merged

Loosen signature in triangular solver from Strided- to AbstractMatrix. #16205

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
b9e1fa7
Loosen signature in triangular solver from Strided- to AbstractMatrix.
andreasnoack May 5, 2016
48719c5
Further improvements
andreasnoack May 7, 2016
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
6 changes: 5 additions & 1 deletion base/linalg/bidiag.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -224,7 +224,6 @@ 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)


BiTriSym = Union{Bidiagonal, Tridiagonal, SymTridiagonal}
BiTri = Union{Bidiagonal, Tridiagonal}
A_mul_B!(C::AbstractMatrix, A::SymTridiagonal, B::BiTriSym) = A_mul_B_td!(C, A, B)
Expand Down Expand Up @@ -361,6 +360,11 @@ function A_mul_B_td!(C::AbstractMatrix, A::AbstractMatrix, B::BiTriSym)
C
end

SpecialMatrix = Union{Bidiagonal, SymTridiagonal, Tridiagonal}
# to avoid ambiguity warning, but shouldn't be necessary
*(A::AbstractTriangular, B::SpecialMatrix) = full(A) * full(B)
*(A::SpecialMatrix, B::SpecialMatrix) = full(A) * full(B)

#Generic multiplication
for func in (:*, :Ac_mul_B, :A_mul_Bc, :/, :A_rdiv_Bc)
@eval ($func){T}(A::Bidiagonal{T}, B::AbstractVector{T}) = ($func)(full(A), B)
Expand Down
21 changes: 19 additions & 2 deletions base/linalg/diagonal.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -104,9 +104,26 @@ end
/{T<:Number}(D::Diagonal, x::T) = Diagonal(D.diag / x)
*(Da::Diagonal, Db::Diagonal) = Diagonal(Da.diag .* Db.diag)
*(D::Diagonal, V::AbstractVector) = D.diag .* V
*(A::AbstractMatrix, D::Diagonal) =
# To avoid ambiguity in the definitions below
for uplo in (:LowerTriangular, :UpperTriangular)
@eval begin
(*)(A::$uplo, D::Diagonal) = $uplo(A.data * D)

function (*)(A::$(Symbol(:Unit, uplo)), D::Diagonal)
B = A.data * D
for i = 1:size(A, 1)
B[i,i] = D.diag[i]
end
return $uplo(B)
end
end
end
(*)(A::AbstractTriangular, D::Diagonal) = error("this method should never be reached")
(*)(D::Diagonal, A::AbstractTriangular) = error("this method should never be reached")
Copy link
Contributor

Choose a reason for hiding this comment

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

are these only necessary to resolve ambiguity?

Copy link
Member Author

Choose a reason for hiding this comment

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

Yes.


(*)(A::AbstractMatrix, D::Diagonal) =
scale!(similar(A, promote_op(*, eltype(A), eltype(D.diag))), A, D.diag)
*(D::Diagonal, A::AbstractMatrix) =
(*)(D::Diagonal, A::AbstractMatrix) =
scale!(similar(A, promote_op(*, eltype(A), eltype(D.diag))), D.diag, A)

A_mul_B!(A::Diagonal,B::AbstractMatrix) = scale!(A.diag,B)
Expand Down
237 changes: 206 additions & 31 deletions base/linalg/triangular.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1291,87 +1291,262 @@ function A_rdiv_Bt!(A::StridedMatrix, B::UnitLowerTriangular)
A
end

for f in (:A_mul_B!, :A_ldiv_B!)
@eval begin
$f(A::UpperTriangular, B::UpperTriangular) =
UpperTriangular($f(A, triu!(B.data)))
$f(A::UnitUpperTriangular, B::UpperTriangular) =
UpperTriangular($f(A, triu!(B.data)))
$f(A::UpperTriangular, B::UnitUpperTriangular) =
UpperTriangular($f(A, triu!(B.data)))
Copy link
Contributor

@tkelman tkelman May 19, 2016 •
edited
Loading

Choose a reason for hiding this comment

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

I don't think this is correct, it uses the diagonals from B.data without setting them to 1 first.

 _/ |\__'_|_|_|\__'_|  |  Commit 1ce4ef0* (0 days old master)
|__/                   |  x86_64-w64-mingw32

julia> A = UpperTriangular(rand(5,5))
5×5 UpperTriangular{Float64,Array{Float64,2}}:
 0.817395  0.512978  0.151383  0.545986   0.680636
  ⋅        0.371017  0.326811  0.0279858  0.242921
  ⋅         ⋅        0.738791  0.0559057  0.272851
  ⋅         ⋅         ⋅        0.780668   0.0159331
  ⋅         ⋅         ⋅         ⋅         0.285675

julia> B = LinAlg.UnitUpperTriangular(rand(5,5))
5×5 Base.LinAlg.UnitUpperTriangular{Float64,Array{Float64,2}}:
 1.0  0.881164  0.0395179   0.265299  0.901799
 0.0  1.0       0.00318107  0.382759  0.909275
 0.0  0.0       1.0         0.525351  0.858504
 0.0  0.0       0.0         1.0       0.447239
 0.0  0.0       0.0         0.0       1.0

julia> B.data
5×5 Array{Float64,2}:
 0.169828  0.881164   0.0395179   0.265299  0.901799
 0.244341  0.701702   0.00318107  0.382759  0.909275
 0.952765  0.907991   0.68711     0.525351  0.858504
 0.924607  0.821667   0.173782    0.737463  0.447239
 0.22737   0.0277924  0.332617    0.311489  0.0209441

julia> full(A) * full(B)
5×5 Array{Float64,2}:
 0.817395  1.23324   0.185317  1.03872   2.25835
 0.0       0.371017  0.327992  0.341686  0.873363
 0.0       0.0       0.738791  0.44403   0.93211
 0.0       0.0       0.0       0.780668  0.365078
 0.0       0.0       0.0       0.0       0.285675

julia> C = LinAlg.A_mul_B!(A, B)
5×5 UpperTriangular{Float64,Array{Float64,2}}:
 0.138816  1.08022   0.137951  0.895375  1.59197
  ⋅        0.260344  0.225736  0.334339  0.63553
  ⋅         ⋅        0.507631  0.429353  0.664973
  ⋅         ⋅         ⋅        0.575714  0.349478
  ⋅         ⋅         ⋅         ⋅        0.0059832

edited to do full(A) * full(B) before mutating

Copy link
Member Author

Choose a reason for hiding this comment

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

True. Time to add some tests to triangular.jl.

$f(A::UnitUpperTriangular, B::UnitUpperTriangular) =
UnitUpperTriangular($f(A, triu!(B.data)))
$f(A::LowerTriangular, B::LowerTriangular) =
LowerTriangular($f(A, tril!(B.data)))
$f(A::UnitLowerTriangular, B::LowerTriangular) =
LowerTriangular($f(A, tril!(B.data)))
$f(A::LowerTriangular, B::UnitLowerTriangular) =
LowerTriangular($f(A, tril!(B.data)))
$f(A::UnitLowerTriangular, B::UnitLowerTriangular) =
LowerTriangular($f(A, tril!(B.data)))
end
end

for f in (:Ac_mul_B!, :At_mul_B!, :Ac_ldiv_B!, :At_ldiv_B!)
@eval begin
$f(A::Union{LowerTriangular,UnitLowerTriangular}, B::UpperTriangular) =
UpperTriangular($f(A, triu!(B.data)))
$f(A::Union{UpperTriangular,UnitUpperTriangular}, B::LowerTriangular) =
LowerTriangular($f(A, tril!(B.data)))
end
end

A_rdiv_B!(A::UpperTriangular, B::Union{UpperTriangular,UnitUpperTriangular}) =
UpperTriangular(A_rdiv_B!(triu!(A.data), B))
A_rdiv_B!(A::LowerTriangular, B::Union{LowerTriangular,UnitLowerTriangular}) =
LowerTriangular(A_rdiv_B!(tril!(A.data), B))

for f in (:A_mul_Bc!, :A_mul_Bt!, :A_rdiv_Bc!, :A_rdiv_Bt!)
@eval begin
$f(A::UpperTriangular, B::Union{LowerTriangular,UnitLowerTriangular}) =
UpperTriangular($f(triu!(A.data), B))
$f(A::LowerTriangular, B::Union{UpperTriangular,UnitUpperTriangular}) =
LowerTriangular($f(tril!(A.data), B))
end
end

# Promotion
## Promotion methods in matmul don't apply to triangular multiplication since it is inplace. Hence we have to make very similar definitions, but without allocation of a result array. For multiplication and unit diagonal division the element type doesn't have to be stable under division whereas that is necessary in the general triangular solve problem.

## Some Triangular-Triangular cases. We might want to write taylored methods for these cases, but I'm not sure it is worth it.
for t in (UpperTriangular, UnitUpperTriangular, LowerTriangular, UnitLowerTriangular)
@eval begin
*(A::Tridiagonal, B::$t) = A_mul_B!(full(A), B)
(*)(A::Tridiagonal, B::$t) = A_mul_B!(full(A), B)
end
end

for f in (:*, :Ac_mul_B, :At_mul_B, :\, :Ac_ldiv_B, :At_ldiv_B)
for (f1, f2) in ((:*, :A_mul_B!), (:\, :A_ldiv_B))
@eval begin
($f)(A::AbstractTriangular, B::AbstractTriangular) = ($f)(A, full(B))
function $f1(A::LowerTriangular, B::LowerTriangular)
TAB = typeof(zero(eltype(A))*zero(eltype(B)) + zero(eltype(A))*zero(eltype(B)))
return LowerTriangular($f2(convert(AbstractMatrix{TAB}, A), copy_oftype(B, TAB)))
end

function $f1(A::UnitLowerTriangular, B::LowerTriangular)
TAB = typeof(zero(eltype(A))*zero(eltype(B)) + zero(eltype(A))*zero(eltype(B)))
return LowerTriangular($f2(convert(AbstractMatrix{TAB}, A), copy_oftype(B, TAB)))
end

function $f1(A::UpperTriangular, B::UpperTriangular)
TAB = typeof(zero(eltype(A))*zero(eltype(B)) + zero(eltype(A))*zero(eltype(B)))
return UpperTriangular($f2(convert(AbstractMatrix{TAB}, A), copy_oftype(B, TAB)))
end

function $f1(A::UnitUpperTriangular, B::UpperTriangular)
TAB = typeof(zero(eltype(A))*zero(eltype(B)) + zero(eltype(A))*zero(eltype(B)))
return UpperTriangular($f2(convert(AbstractMatrix{TAB}, A), copy_oftype(B, TAB)))
end
end
end
for f in (:A_mul_Bc, :A_mul_Bt, :Ac_mul_Bc, :At_mul_Bt, :/, :A_rdiv_Bc, :A_rdiv_Bt)

for (f1, f2) in ((:Ac_mul_B, :Ac_mul_B!), (:At_mul_B, :At_mul_B!),
(:Ac_ldiv_B, Ac_ldiv_B!), (:At_ldiv_B, :At_ldiv_B!))
@eval begin
($f)(A::AbstractTriangular, B::AbstractTriangular) = ($f)(full(A), B)
function $f1(A::UpperTriangular, B::LowerTriangular)
TAB = typeof(zero(eltype(A))*zero(eltype(B)) + zero(eltype(A))*zero(eltype(B)))
return LowerTriangular($f2(convert(AbstractMatrix{TAB}, A), copy_oftype(B, TAB)))
end

function $f1(A::UnitUpperTriangular, B::LowerTriangular)
TAB = typeof(zero(eltype(A))*zero(eltype(B)) + zero(eltype(A))*zero(eltype(B)))
return LowerTriangular($f2(convert(AbstractMatrix{TAB}, A), copy_oftype(B, TAB)))
end

function $f1(A::LowerTriangular, B::UpperTriangular)
TAB = typeof(zero(eltype(A))*zero(eltype(B)) + zero(eltype(A))*zero(eltype(B)))
return UpperTriangular($f2(convert(AbstractMatrix{TAB}, A), copy_oftype(B, TAB)))
end

function $f1(A::UnitLowerTriangular, B::UpperTriangular)
TAB = typeof(zero(eltype(A))*zero(eltype(B)) + zero(eltype(A))*zero(eltype(B)))
return UpperTriangular($f2(convert(AbstractMatrix{TAB}, A), copy_oftype(B, TAB)))
end
end
end

function (/)(A::LowerTriangular, B::LowerTriangular)
TAB = typeof(zero(eltype(A))*zero(eltype(B)) + zero(eltype(A))*zero(eltype(B))/
one(eltype(A)))
return LowerTriangular(A_rdiv_B!(copy_oftype(A, TAB), convert(AbstractMatrix{TAB}, B)))
end
function (/)(A::LowerTriangular, B::UnitLowerTriangular)
TAB = typeof(zero(eltype(A))*zero(eltype(B)) + zero(eltype(A))*zero(eltype(B)))
return LowerTriangular(A_rdiv_B!(copy_oftype(A, TAB), convert(AbstractMatrix{TAB}, B)))
end
function (/)(A::UpperTriangular, B::UpperTriangular)
TAB = typeof(zero(eltype(A))*zero(eltype(B)) + zero(eltype(A))*zero(eltype(B))/
one(eltype(A)))
return UpperTriangular(A_rdiv_B!(copy_oftype(A, TAB), convert(AbstractMatrix{TAB}, B)))
end
function (/)(A::UpperTriangular, B::UnitUpperTriangular)
TAB = typeof(zero(eltype(A))*zero(eltype(B)) + zero(eltype(A))*zero(eltype(B)))
return UpperTriangular(A_rdiv_B!(copy_oftype(A, TAB), convert(AbstractMatrix{TAB}, B)))
end

for (f1, f2) in ((:A_mul_Bc, :A_mul_Bc!), (:A_mul_Bt, :A_mul_Bt!),
(:A_rdiv_Bc, :A_rdiv_Bc!), (:A_rdiv_Bt, :A_rdiv_Bt!))
@eval begin
function $f1(A::LowerTriangular, B::UpperTriangular)
TAB = typeof(zero(eltype(A))*zero(eltype(B)) + zero(eltype(A))*zero(eltype(B)))
return LowerTriangular($f2(copy_oftype(A, TAB), convert(AbstractMatrix{TAB}, B)))
end

function $f1(A::LowerTriangular, B::UnitUpperTriangular)
TAB = typeof(zero(eltype(A))*zero(eltype(B)) + zero(eltype(A))*zero(eltype(B)))
return LowerTriangular($f2(copy_oftype(A, TAB), convert(AbstractMatrix{TAB}, B)))
end

function $f1(A::UpperTriangular, B::LowerTriangular)
TAB = typeof(zero(eltype(A))*zero(eltype(B)) + zero(eltype(A))*zero(eltype(B)))
return UpperTriangular($f2(copy_oftype(A, TAB), convert(AbstractMatrix{TAB}, B)))
end

function $f1(A::UpperTriangular, B::UnitLowerTriangular)
TAB = typeof(zero(eltype(A))*zero(eltype(B)) + zero(eltype(A))*zero(eltype(B)))
return UpperTriangular($f2(copy_oftype(A, TAB), convert(AbstractMatrix{TAB}, B)))
end
end
end

## The general promotion methods

for (f, g) in ((:*, :A_mul_B!), (:Ac_mul_B, :Ac_mul_B!), (:At_mul_B, :At_mul_B!))
@eval begin
function ($f)(A::AbstractTriangular, B::AbstractTriangular)
TAB = typeof(zero(eltype(A))*zero(eltype(B)) + zero(eltype(A))*zero(eltype(B)))
BB = similar(B, TAB, size(B))
copy!(BB, B)
($g)(convert(AbstractArray{TAB}, A), BB)
end
end
end
for (f, g) in ((:A_mul_Bc, :A_mul_Bc!), (:A_mul_Bt, :A_mul_Bt!))
@eval begin
function ($f)(A::AbstractTriangular, B::AbstractTriangular)
TAB = typeof(zero(eltype(A))*zero(eltype(B)) + zero(eltype(A))*zero(eltype(B)))
AA = similar(A, TAB, size(A))
copy!(AA, A)
($g)(AA, convert(AbstractArray{TAB}, B))
end
end
end

for mat in (:AbstractVector, :AbstractMatrix)

### 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::AbstractTriangular{TA}, B::StridedVecOrMat{TB})
TAB = typeof(zero(TA)*zero(TB) + zero(TA)*zero(TB))
($g)(convert(AbstractArray{TAB}, A), copy_oftype(B, TAB))
function ($f)(A::AbstractTriangular, B::$mat)
TAB = typeof(zero(eltype(A))*zero(eltype(B)) + zero(eltype(A))*zero(eltype(B)))
BB = similar(B, TAB, size(B))
copy!(BB, B)
($g)(convert(AbstractArray{TAB}, A), BB)
end
end
end
### Left division with triangle to the left hence rhs cannot be transposed. No quotients.
for (f, g) in ((:\, :A_ldiv_B!), (:Ac_ldiv_B, :Ac_ldiv_B!), (:At_ldiv_B, :At_ldiv_B!))
@eval begin
function ($f){TA,TB,S}(A::Union{UnitUpperTriangular{TA,S},UnitLowerTriangular{TA,S}}, B::StridedVecOrMat{TB})
TAB = typeof(zero(TA)*zero(TB) + zero(TA)*zero(TB))
($g)(convert(AbstractArray{TAB}, A), copy_oftype(B, TAB))
function ($f)(A::Union{UnitUpperTriangular,UnitLowerTriangular}, B::$mat)
TAB = typeof(zero(eltype(A))*zero(eltype(B)) + zero(eltype(A))*zero(eltype(B)))
BB = similar(B, TAB, size(B))
copy!(BB, B)
($g)(convert(AbstractArray{TAB}, A), BB)
end
end
end
### Left division with triangle to the left hence rhs cannot be transposed. Quotients.
for (f, g) in ((:\, :A_ldiv_B!), (:Ac_ldiv_B, :Ac_ldiv_B!), (:At_ldiv_B, :At_ldiv_B!))
@eval begin
function ($f){TA,TB,S}(A::Union{UpperTriangular{TA,S},LowerTriangular{TA,S}}, B::StridedVecOrMat{TB})
TAB = typeof((zero(TA)*zero(TB) + zero(TA)*zero(TB))/one(TA))
($g)(convert(AbstractArray{TAB}, A), copy_oftype(B, TAB))
function ($f)(A::Union{UpperTriangular,LowerTriangular}, B::$mat)
TAB = typeof((zero(eltype(A))*zero(eltype(B)) + zero(eltype(A))*zero(eltype(B)))/one(eltype(A)))
BB = similar(B, TAB, size(B))
copy!(BB, B)
($g)(convert(AbstractArray{TAB}, A), BB)
end
end
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::AbstractTriangular{TB})
TAB = typeof(zero(TA)*zero(TB) + zero(TA)*zero(TB))
($g)(copy_oftype(A, TAB), convert(AbstractArray{TAB}, B))
function ($f)(A::$mat, B::AbstractTriangular)
TAB = typeof(zero(eltype(A))*zero(eltype(B)) + zero(eltype(A))*zero(eltype(B)))
AA = similar(A, TAB, size(A))
copy!(AA, A)
($g)(AA, convert(AbstractArray{TAB}, B))
end
end
end
### Right division with triangle to the right hence lhs cannot be transposed. No quotients.
for (f, g) in ((:/, :A_rdiv_B!), (:A_rdiv_Bc, :A_rdiv_Bc!), (:A_rdiv_Bt, :A_rdiv_Bt!))
@eval begin
function ($f){TA,TB,S}(A::StridedVecOrMat{TA}, B::Union{UnitUpperTriangular{TB,S},UnitLowerTriangular{TB,S}})
TAB = typeof(zero(TA)*zero(TB) + zero(TA)*zero(TB))
($g)(copy_oftype(A, TAB), convert(AbstractArray{TAB}, B))
function ($f)(A::$mat, B::Tuple{UnitUpperTriangular, UnitLowerTriangular})
TAB = typeof(zero(eltype(A))*zero(eltype(B)) + zero(eltype(A))*zero(eltype(B)))
AA = similar(A, TAB, size(A))
copy!(AA, A)
($g)(AA, convert(AbstractArray{TAB}, B))
end
end
end

### Right division with triangle to the right hence lhs cannot be transposed. Quotients.
for (f, g) in ((:/, :A_rdiv_B!), (:A_rdiv_Bc, :A_rdiv_Bc!), (:A_rdiv_Bt, :A_rdiv_Bt!))
@eval begin
function ($f){TA,TB,S}(A::StridedVecOrMat{TA}, B::Union{UpperTriangular{TB,S},LowerTriangular{TB,S}})
TAB = typeof((zero(TA)*zero(TB) + zero(TA)*zero(TB))/one(TA))
($g)(copy_oftype(A, TAB), convert(AbstractArray{TAB}, B))
end
end
end
### Fallbacks brought in from linalg/bidiag.jl while fixing #14506.
# Eventually the above promotion methods should be generalized as
# was done for bidiagonal matrices in #14506.
At_ldiv_B(A::AbstractTriangular, B::AbstractVecOrMat) = At_ldiv_B!(A, copy(B))
Ac_ldiv_B(A::AbstractTriangular, B::AbstractVecOrMat) = Ac_ldiv_B!(A, copy(B))
function ($f)(A::$mat, B::Union{UpperTriangular,LowerTriangular})
TAB = typeof((zero(eltype(A))*zero(eltype(B)) + zero(eltype(A))*zero(eltype(B)))/one(eltype(A)))
AA = similar(A, TAB, size(A))
copy!(AA, A)
($g)(AA, convert(AbstractArray{TAB}, B))
end
end
end
end

# If these are not defined, the they will fallback to the versions in matmul.jl
# and dispatch to generic_matmatmul! which is very costly to compile. The methods
# below might compute an unnecessary copy. Eliminating the copy requires adding
# all the promotion logic here once again. Since these methods are probably relatively
# rare, we chose not to bother for now.
Ac_mul_B(A::AbstractMatrix, B::AbstractTriangular) = (*)(ctranspose(A), B)
At_mul_B(A::AbstractMatrix, B::AbstractTriangular) = (*)(transpose(A), B)
A_mul_Bc(A::AbstractTriangular, B::AbstractMatrix) = (*)(A, ctranspose(B))
A_mul_Bt(A::AbstractTriangular, B::AbstractMatrix) = (*)(A, transpose(B))
Ac_mul_Bc(A::AbstractTriangular, B::AbstractTriangular) = Ac_mul_B(A, B')
Ac_mul_Bc(A::AbstractTriangular, B::AbstractMatrix) = Ac_mul_B(A, B')
Ac_mul_Bc(A::AbstractMatrix, B::AbstractTriangular) = A_mul_Bc(A', B)
At_mul_Bt(A::AbstractTriangular, B::AbstractTriangular) = At_mul_B(A, B.')
At_mul_Bt(A::AbstractTriangular, B::AbstractMatrix) = At_mul_B(A, B.')
At_mul_Bt(A::AbstractMatrix, B::AbstractTriangular) = A_mul_Bt(A.', B)

# Complex matrix logarithm for the upper triangular factor, see:
# Al-Mohy and Higham, "Improved inverse scaling and squaring algorithms for
Expand Down
4 changes: 2 additions & 2 deletions base/sparse/sparsevector.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1525,7 +1525,7 @@ for isunittri in (true, false), islowertri in (true, false)
(true, :(Ac_ldiv_B), :(Ac_ldiv_B!)) )

# broad method where elements are Numbers
@eval function ($func){TA<:Number,Tb<:Number,S}(A::$tritype{TA,S}, b::SparseVector{Tb})
@eval function ($func){TA<:Number,Tb<:Number,S<:AbstractMatrix}(A::$tritype{TA,S}, b::SparseVector{Tb})
TAb = $(isunittri ?
:(typeof(zero(TA)*zero(Tb) + zero(TA)*zero(Tb))) :
:(typeof((zero(TA)*zero(Tb) + zero(TA)*zero(Tb))/one(TA))) )
Expand All @@ -1551,7 +1551,7 @@ for isunittri in (true, false), islowertri in (true, false)
end

# fallback where elements are not Numbers
@eval ($func){TA,Tb,S}(A::$tritype{TA,S}, b::SparseVector{Tb}) = ($ipfunc)(A, copy(b))
@eval ($func)(A::$tritype, b::SparseVector) = ($ipfunc)(A, copy(b))
end

# build in-place left-division operations
Expand Down
4 changes: 4 additions & 0 deletions test/linalg/triangular.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -24,6 +24,7 @@ for elty1 in (Float32, Float64, BigFloat, Complex64, Complex128, Complex{BigFloa
# Construct test matrix
A1 = t1(elty1 == Int ? rand(1:7, n, n) : convert(Matrix{elty1}, (elty1 <: Complex ? complex(randn(n, n), randn(n, n)) : randn(n, n)) |> t -> chol(t't) |> t -> uplo1 == :U ? t : ctranspose(t)))


debug && println("elty1: $elty1, A1: $t1")

# Convert
Expand Down Expand Up @@ -489,3 +490,6 @@ end
# Test that UpperTriangular(LowerTriangular) throws. See #16201
@test_throws ArgumentError LowerTriangular(UpperTriangular(randn(3,3)))
@test_throws ArgumentError UpperTriangular(LowerTriangular(randn(3,3)))

# Issue 16196
@test UpperTriangular(eye(3)) \ sub(ones(3), [1,2,3]) == ones(3)