Skip to content

Commit

Permalink
Linalg: Reduce allocations in triangular tests (#53634)
Browse files Browse the repository at this point in the history 
Reuses a pre-allocated matrix in tests to avoid allocating a fresh
matrix in every call.
  • Loading branch information
@jishnub
jishnub authored Mar 9, 2024
1 parent bb35dc9 commit 53048b2
Showing 1 changed file with 108 additions and 107 deletions.
215 changes: 108 additions & 107 deletions stdlib/LinearAlgebra/test/triangular.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -35,6 +35,7 @@ debug && println("Test basic type functionality")

# 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 -> cholesky(t't).U |> t -> uplo1 === :U ? t : copy(t')))
M1 = Matrix(A1)
@test t1(A1) === A1
@test t1{elty1}(A1) === A1
# test the ctor works for AbstractMatrix
Expand Down Expand Up @@ -68,7 +69,7 @@ debug && println("Test basic type functionality")
@test simA1 == A1

# getindex
let mA1 = Matrix(A1)
let mA1 = M1
# linear indexing
for i in 1:length(A1)
@test A1[i] == mA1[i]
Expand Down Expand Up @@ -141,8 +142,8 @@ debug && println("Test basic type functionality")
#tril/triu
if uplo1 === :L
@test tril(A1,0) == A1
@test tril(A1,-1) == LowerTriangular(tril(Matrix(A1), -1))
@test tril(A1,1) == t1(tril(tril(Matrix(A1), 1)))
@test tril(A1,-1) == LowerTriangular(tril(M1, -1))
@test tril(A1,1) == t1(tril(tril(M1, 1)))
@test tril(A1, -n - 2) == zeros(size(A1))
@test tril(A1, n) == A1
@test triu(A1,0) == t1(diagm(0 => diag(A1)))
Expand All @@ -152,8 +153,8 @@ debug && println("Test basic type functionality")
@test triu(A1, n + 2) == zeros(size(A1))
else
@test triu(A1,0) == A1
@test triu(A1,1) == UpperTriangular(triu(Matrix(A1), 1))
@test triu(A1,-1) == t1(triu(triu(Matrix(A1), -1)))
@test triu(A1,1) == UpperTriangular(triu(M1, 1))
@test triu(A1,-1) == t1(triu(triu(M1, -1)))
@test triu(A1, -n) == A1
@test triu(A1, n + 2) == zeros(size(A1))
@test tril(A1,0) == t1(diagm(0 => diag(A1)))
Expand All @@ -169,10 +170,10 @@ debug && println("Test basic type functionality")
# [c]transpose[!] (test views as well, see issue #14317)
let vrange = 1:n-1, viewA1 = t1(view(A1.data, vrange, vrange))
# transpose
@test copy(transpose(A1)) == transpose(Matrix(A1))
@test copy(transpose(A1)) == transpose(M1)
@test copy(transpose(viewA1)) == transpose(Matrix(viewA1))
# adjoint
@test copy(A1') == Matrix(A1)'
@test copy(A1') == M1'
@test copy(viewA1') == Matrix(viewA1)'
# transpose!
@test transpose!(copy(A1)) == transpose(A1)
Expand All @@ -185,28 +186,28 @@ debug && println("Test basic type functionality")
end

# diag
@test diag(A1) == diag(Matrix(A1))
@test diag(A1) == diag(M1)

# tr
@test tr(A1)::elty1 == tr(Matrix(A1))
@test tr(A1)::elty1 == tr(M1)

# real
@test real(A1) == real(Matrix(A1))
@test imag(A1) == imag(Matrix(A1))
@test abs.(A1) == abs.(Matrix(A1))
@test real(A1) == real(M1)
@test imag(A1) == imag(M1)
@test abs.(A1) == abs.(M1)

# zero
if A1 isa UpperTriangular || A1 isa LowerTriangular
@test zero(A1) == zero(parent(A1))
end

# Unary operations
@test -A1 == -Matrix(A1)
@test -A1 == -M1

# copy and copyto! (test views as well, see issue #14317)
let vrange = 1:n-1, viewA1 = t1(view(A1.data, vrange, vrange))
# copy
@test copy(A1) == copy(Matrix(A1))
@test copy(A1) == copy(M1)
@test copy(viewA1) == copy(Matrix(viewA1))
# copyto!
B = similar(A1)
Expand Down Expand Up @@ -283,25 +284,25 @@ debug && println("Test basic type functionality")
end

# Binary operations
@test A1*0.5 == Matrix(A1)*0.5
@test 0.5*A1 == 0.5*Matrix(A1)
@test A1/0.5 == Matrix(A1)/0.5
@test 0.5\A1 == 0.5\Matrix(A1)
@test A1*0.5 == M1*0.5
@test 0.5*A1 == 0.5*M1
@test A1/0.5 == M1/0.5
@test 0.5\A1 == 0.5\M1

# inversion
@test inv(A1) ≈ inv(lu(Matrix(A1)))
inv(Matrix(A1)) # issue #11298
@test inv(A1) ≈ inv(lu(M1))
inv(M1) # issue #11298
@test isa(inv(A1), t1)
# make sure the call to LAPACK works right
if elty1 <: BlasFloat
@test LinearAlgebra.inv!(copy(A1)) ≈ inv(lu(Matrix(A1)))
@test LinearAlgebra.inv!(copy(A1)) ≈ inv(lu(M1))
end

# Determinant
@test det(A1) ≈ det(lu(Matrix(A1))) atol=sqrt(eps(real(float(one(elty1)))))*n*n
@test logdet(A1) ≈ logdet(lu(Matrix(A1))) atol=sqrt(eps(real(float(one(elty1)))))*n*n
@test det(A1) ≈ det(lu(M1)) atol=sqrt(eps(real(float(one(elty1)))))*n*n
@test logdet(A1) ≈ logdet(lu(M1)) atol=sqrt(eps(real(float(one(elty1)))))*n*n
lada, ladb = logabsdet(A1)
flada, fladb = logabsdet(lu(Matrix(A1)))
flada, fladb = logabsdet(lu(M1))
@test lada ≈ flada atol=sqrt(eps(real(float(one(elty1)))))*n*n
@test ladb ≈ fladb atol=sqrt(eps(real(float(one(elty1)))))*n*n

Expand All @@ -324,7 +325,7 @@ debug && println("Test basic type functionality")
for p in (1.0, Inf)
@test cond(A1,p) ≈ cond(A1,p) atol=(cond(A1,p)+cond(A1,p))
end
@test cond(A1,2) == cond(Matrix(A1),2)
@test cond(A1,2) == cond(M1,2)
end

if !(elty1 in (BigFloat, Complex{BigFloat})) # Not implemented yet
Expand All @@ -333,9 +334,9 @@ debug && println("Test basic type functionality")
svdvals(A1)
end

@test ((A1*A1)::t1) ≈ Matrix(A1) * Matrix(A1)
@test ((A1/A1)::t1) ≈ Matrix(A1) / Matrix(A1)
@test ((A1\A1)::t1) ≈ Matrix(A1) \ Matrix(A1)
@test ((A1*A1)::t1) ≈ M1 * M1
@test ((A1/A1)::t1) ≈ M1 / M1
@test ((A1\A1)::t1) ≈ M1 \ M1

# Begin loop for second Triangular matrix
for elty2 in (Float32, Float64, BigFloat, ComplexF32, ComplexF64, Complex{BigFloat}, Int)
Expand All @@ -347,7 +348,7 @@ debug && println("Test basic type functionality")
debug && println("elty1: $elty1, A1: $t1, elty2: $elty2, A2: $t2")

A2 = t2(elty2 == Int ? rand(1:7, n, n) : convert(Matrix{elty2}, (elty2 <: Complex ? complex.(randn(n, n), randn(n, n)) : randn(n, n)) |> t -> cholesky(t't).U |> t -> uplo2 === :U ? t : copy(t')))

M2 = Matrix(A2)
# Convert
if elty1 <: Real && !(elty2 <: Integer)
@test convert(AbstractMatrix{elty2}, A1) == t1(convert(Matrix{elty2}, A1.data))
Expand All @@ -356,21 +357,21 @@ debug && println("Test basic type functionality")
end

# Binary operations
@test A1 + A2 == Matrix(A1) + Matrix(A2)
@test A1 - A2 == Matrix(A1) - Matrix(A2)
@test A1 + A2 == M1 + M2
@test A1 - A2 == M1 - M2

# Triangular-Triangular multiplication and division
@test A1*A2 ≈ Matrix(A1)*Matrix(A2)
@test transpose(A1)*A2 ≈ transpose(Matrix(A1))*Matrix(A2)
@test transpose(A1)*adjoint(A2) ≈ transpose(Matrix(A1))*adjoint(Matrix(A2))
@test adjoint(A1)*transpose(A2) ≈ adjoint(Matrix(A1))*transpose(Matrix(A2))
@test A1'A2 ≈ Matrix(A1)'Matrix(A2)
@test A1*transpose(A2) ≈ Matrix(A1)*transpose(Matrix(A2))
@test A1*A2' ≈ Matrix(A1)*Matrix(A2)'
@test transpose(A1)*transpose(A2) ≈ transpose(Matrix(A1))*transpose(Matrix(A2))
@test A1'A2' ≈ Matrix(A1)'Matrix(A2)'
@test A1/A2 ≈ Matrix(A1)/Matrix(A2)
@test A1\A2 ≈ Matrix(A1)\Matrix(A2)
@test A1*A2 ≈ M1*M2
@test transpose(A1)*A2 ≈ transpose(M1)*M2
@test transpose(A1)*adjoint(A2) ≈ transpose(M1)*adjoint(M2)
@test adjoint(A1)*transpose(A2) ≈ adjoint(M1)*transpose(M2)
@test A1'A2 ≈ M1'M2
@test A1*transpose(A2) ≈ M1*transpose(M2)
@test A1*A2' ≈ M1*M2'
@test transpose(A1)*transpose(A2) ≈ transpose(M1)*transpose(M2)
@test A1'A2' ≈ M1'M2'
@test A1/A2 ≈ M1/M2
@test A1\A2 ≈ M1\M2
if uplo1 === :U && uplo2 === :U
if t1 === UnitUpperTriangular && t2 === UnitUpperTriangular
@test A1*A2 isa UnitUpperTriangular
Expand Down Expand Up @@ -411,20 +412,20 @@ debug && println("Test basic type functionality")
@test_throws DimensionMismatch A2' * offsizeA
@test_throws DimensionMismatch A2 * offsizeA
if (uplo1 == uplo2 && elty1 == elty2 != Int && t1 != UnitLowerTriangular && t1 != UnitUpperTriangular)
@test rdiv!(copy(A1), A2)::t1 ≈ A1/A2 ≈ Matrix(A1)/Matrix(A2)
@test ldiv!(A2, copy(A1))::t1 ≈ A2\A1 ≈ Matrix(A2)\Matrix(A1)
@test rdiv!(copy(A1), A2)::t1 ≈ A1/A2 ≈ M1/M2
@test ldiv!(A2, copy(A1))::t1 ≈ A2\A1 ≈ M2\M1
end
if (uplo1 != uplo2 && elty1 == elty2 != Int && t2 != UnitLowerTriangular && t2 != UnitUpperTriangular)
@test lmul!(adjoint(A1), copy(A2)) ≈ A1'*A2 ≈ Matrix(A1)'*Matrix(A2)
@test lmul!(transpose(A1), copy(A2)) ≈ transpose(A1)*A2 ≈ transpose(Matrix(A1))*Matrix(A2)
@test ldiv!(adjoint(A1), copy(A2)) ≈ A1'\A2 ≈ Matrix(A1)'\Matrix(A2)
@test ldiv!(transpose(A1), copy(A2)) ≈ transpose(A1)\A2 ≈ transpose(Matrix(A1))\Matrix(A2)
@test lmul!(adjoint(A1), copy(A2)) ≈ A1'*A2 ≈ M1'*M2
@test lmul!(transpose(A1), copy(A2)) ≈ transpose(A1)*A2 ≈ transpose(M1)*M2
@test ldiv!(adjoint(A1), copy(A2)) ≈ A1'\A2 ≈ M1'\M2
@test ldiv!(transpose(A1), copy(A2)) ≈ transpose(A1)\A2 ≈ transpose(M1)\M2
end
if (uplo1 != uplo2 && elty1 == elty2 != Int && t1 != UnitLowerTriangular && t1 != UnitUpperTriangular)
@test rmul!(copy(A1), adjoint(A2)) ≈ A1*A2' ≈ Matrix(A1)*Matrix(A2)'
@test rmul!(copy(A1), transpose(A2)) ≈ A1*transpose(A2) ≈ Matrix(A1)*transpose(Matrix(A2))
@test rdiv!(copy(A1), adjoint(A2)) ≈ A1/A2' ≈ Matrix(A1)/Matrix(A2)'
@test rdiv!(copy(A1), transpose(A2)) ≈ A1/transpose(A2) ≈ Matrix(A1)/transpose(Matrix(A2))
@test rmul!(copy(A1), adjoint(A2)) ≈ A1*A2' ≈ M1*M2'
@test rmul!(copy(A1), transpose(A2)) ≈ A1*transpose(A2) ≈ M1*transpose(M2)
@test rdiv!(copy(A1), adjoint(A2)) ≈ A1/A2' ≈ M1/M2'
@test rdiv!(copy(A1), transpose(A2)) ≈ A1/transpose(A2) ≈ M1/transpose(M2)
end
end
end
Expand All @@ -435,55 +436,55 @@ debug && println("Test basic type functionality")
debug && println("elty1: $elty1, A1: $t1, B: $eltyB")

Tri = Tridiagonal(rand(eltyB,n-1),rand(eltyB,n),rand(eltyB,n-1))
@test lmul!(Tri,copy(A1)) ≈ Tri*Matrix(A1)
@test lmul!(Tri,copy(A1)) ≈ Tri*M1
Tri = Tridiagonal(rand(eltyB,n-1),rand(eltyB,n),rand(eltyB,n-1))
C = Matrix{promote_type(elty1,eltyB)}(undef, n, n)
mul!(C, Tri, A1)
@test C ≈ Tri*Matrix(A1)
@test C ≈ Tri*M1
Tri = Tridiagonal(rand(eltyB,n-1),rand(eltyB,n),rand(eltyB,n-1))
mul!(C, A1, Tri)
@test C ≈ Matrix(A1)*Tri
@test C ≈ M1*Tri

# Triangular-dense Matrix/vector multiplication
@test A1*B[:,1] ≈ Matrix(A1)*B[:,1]
@test A1*B ≈ Matrix(A1)*B
@test transpose(A1)*B[:,1] ≈ transpose(Matrix(A1))*B[:,1]
@test A1'B[:,1] ≈ Matrix(A1)'B[:,1]
@test transpose(A1)*B ≈ transpose(Matrix(A1))*B
@test A1'B ≈ Matrix(A1)'B
@test A1*transpose(B) ≈ Matrix(A1)*transpose(B)
@test adjoint(A1)*transpose(B) ≈ Matrix(A1)'*transpose(B)
@test transpose(A1)*adjoint(B) ≈ transpose(Matrix(A1))*adjoint(B)
@test A1*B' ≈ Matrix(A1)*B'
@test B*A1 ≈ B*Matrix(A1)
@test transpose(B[:,1])*A1 ≈ transpose(B[:,1])*Matrix(A1)
@test B[:,1]'A1 ≈ B[:,1]'Matrix(A1)
@test transpose(B)*A1 ≈ transpose(B)*Matrix(A1)
@test transpose(B)*adjoint(A1) ≈ transpose(B)*Matrix(A1)'
@test adjoint(B)*transpose(A1) ≈ adjoint(B)*transpose(Matrix(A1))
@test B'A1 ≈ B'Matrix(A1)
@test B*transpose(A1) ≈ B*transpose(Matrix(A1))
@test B*A1' ≈ B*Matrix(A1)'
@test transpose(B[:,1])*transpose(A1) ≈ transpose(B[:,1])*transpose(Matrix(A1))
@test B[:,1]'A1' ≈ B[:,1]'Matrix(A1)'
@test transpose(B)*transpose(A1) ≈ transpose(B)*transpose(Matrix(A1))
@test B'A1' ≈ B'Matrix(A1)'
@test A1*B[:,1] ≈ M1*B[:,1]
@test A1*B ≈ M1*B
@test transpose(A1)*B[:,1] ≈ transpose(M1)*B[:,1]
@test A1'B[:,1] ≈ M1'B[:,1]
@test transpose(A1)*B ≈ transpose(M1)*B
@test A1'B ≈ M1'B
@test A1*transpose(B) ≈ M1*transpose(B)
@test adjoint(A1)*transpose(B) ≈ M1'*transpose(B)
@test transpose(A1)*adjoint(B) ≈ transpose(M1)*adjoint(B)
@test A1*B' ≈ M1*B'
@test B*A1 ≈ B*M1
@test transpose(B[:,1])*A1 ≈ transpose(B[:,1])*M1
@test B[:,1]'A1 ≈ B[:,1]'M1
@test transpose(B)*A1 ≈ transpose(B)*M1
@test transpose(B)*adjoint(A1) ≈ transpose(B)*M1'
@test adjoint(B)*transpose(A1) ≈ adjoint(B)*transpose(M1)
@test B'A1 ≈ B'M1
@test B*transpose(A1) ≈ B*transpose(M1)
@test B*A1' ≈ B*M1'
@test transpose(B[:,1])*transpose(A1) ≈ transpose(B[:,1])*transpose(M1)
@test B[:,1]'A1' ≈ B[:,1]'M1'
@test transpose(B)*transpose(A1) ≈ transpose(B)*transpose(M1)
@test B'A1' ≈ B'M1'

if eltyB == elty1
@test mul!(similar(B), A1, B) ≈ Matrix(A1)*B
@test mul!(similar(B), A1, adjoint(B)) ≈ Matrix(A1)*B'
@test mul!(similar(B), A1, transpose(B)) ≈ Matrix(A1)*transpose(B)
@test mul!(similar(B), adjoint(A1), adjoint(B)) ≈ Matrix(A1)'*B'
@test mul!(similar(B), transpose(A1), transpose(B)) ≈ transpose(Matrix(A1))*transpose(B)
@test mul!(similar(B), transpose(A1), adjoint(B)) ≈ transpose(Matrix(A1))*B'
@test mul!(similar(B), adjoint(A1), transpose(B)) ≈ Matrix(A1)'*transpose(B)
@test mul!(similar(B), adjoint(A1), B) ≈ Matrix(A1)'*B
@test mul!(similar(B), transpose(A1), B) ≈ transpose(Matrix(A1))*B
@test mul!(similar(B), A1, B) ≈ M1*B
@test mul!(similar(B), A1, adjoint(B)) ≈ M1*B'
@test mul!(similar(B), A1, transpose(B)) ≈ M1*transpose(B)
@test mul!(similar(B), adjoint(A1), adjoint(B)) ≈ M1'*B'
@test mul!(similar(B), transpose(A1), transpose(B)) ≈ transpose(M1)*transpose(B)
@test mul!(similar(B), transpose(A1), adjoint(B)) ≈ transpose(M1)*B'
@test mul!(similar(B), adjoint(A1), transpose(B)) ≈ M1'*transpose(B)
@test mul!(similar(B), adjoint(A1), B) ≈ M1'*B
@test mul!(similar(B), transpose(A1), B) ≈ transpose(M1)*B
# test also vector methods
B1 = vec(B[1,:])
@test mul!(similar(B1), A1, B1) ≈ Matrix(A1)*B1
@test mul!(similar(B1), adjoint(A1), B1) ≈ Matrix(A1)'*B1
@test mul!(similar(B1), transpose(A1), B1) ≈ transpose(Matrix(A1))*B1
@test mul!(similar(B1), A1, B1) ≈ M1*B1
@test mul!(similar(B1), adjoint(A1), B1) ≈ M1'*B1
@test mul!(similar(B1), transpose(A1), B1) ≈ transpose(M1)*B1
end
#error handling
Ann, Bmm, bm = A1, Matrix{eltyB}(undef, n+1, n+1), Vector{eltyB}(undef, n+1)
Expand All @@ -495,30 +496,30 @@ debug && println("Test basic type functionality")
@test_throws DimensionMismatch rmul!(Bmm, transpose(Ann))

# ... and division
@test A1\B[:,1] ≈ Matrix(A1)\B[:,1]
@test A1\B ≈ Matrix(A1)\B
@test transpose(A1)\B[:,1] ≈ transpose(Matrix(A1))\B[:,1]
@test A1'\B[:,1] ≈ Matrix(A1)'\B[:,1]
@test transpose(A1)\B ≈ transpose(Matrix(A1))\B
@test A1'\B ≈ Matrix(A1)'\B
@test A1\transpose(B) ≈ Matrix(A1)\transpose(B)
@test A1\B' ≈ Matrix(A1)\B'
@test transpose(A1)\transpose(B) ≈ transpose(Matrix(A1))\transpose(B)
@test A1'\B' ≈ Matrix(A1)'\B'
@test A1\B[:,1] ≈ M1\B[:,1]
@test A1\B ≈ M1\B
@test transpose(A1)\B[:,1] ≈ transpose(M1)\B[:,1]
@test A1'\B[:,1] ≈ M1'\B[:,1]
@test transpose(A1)\B ≈ transpose(M1)\B
@test A1'\B ≈ M1'\B
@test A1\transpose(B) ≈ M1\transpose(B)
@test A1\B' ≈ M1\B'
@test transpose(A1)\transpose(B) ≈ transpose(M1)\transpose(B)
@test A1'\B' ≈ M1'\B'
Ann, bm = A1, Vector{elty1}(undef,n+1)
@test_throws DimensionMismatch Ann\bm
@test_throws DimensionMismatch Ann'\bm
@test_throws DimensionMismatch transpose(Ann)\bm
if t1 == UpperTriangular || t1 == LowerTriangular
@test_throws SingularException ldiv!(t1(zeros(elty1, n, n)), fill(eltyB(1), n))
end
@test B/A1 ≈ B/Matrix(A1)
@test B/transpose(A1) ≈ B/transpose(Matrix(A1))
@test B/A1' ≈ B/Matrix(A1)'
@test transpose(B)/A1 ≈ transpose(B)/Matrix(A1)
@test B'/A1 ≈ B'/Matrix(A1)
@test transpose(B)/transpose(A1) ≈ transpose(B)/transpose(Matrix(A1))
@test B'/A1' ≈ B'/Matrix(A1)'
@test B/A1 ≈ B/M1
@test B/transpose(A1) ≈ B/transpose(M1)
@test B/A1' ≈ B/M1'
@test transpose(B)/A1 ≈ transpose(B)/M1
@test B'/A1 ≈ B'/M1
@test transpose(B)/transpose(A1) ≈ transpose(B)/transpose(M1)
@test B'/A1' ≈ B'/M1'

# Error bounds
!(elty1 in (BigFloat, Complex{BigFloat})) && !(eltyB in (BigFloat, Complex{BigFloat})) && errorbounds(A1, A1\B, B)
Expand Down

0 comments on commit 53048b2

Please sign in to comment.