Skip to content

Commit 53048b2

Browse files
authored
Linalg: Reduce allocations in triangular tests (JuliaLang#53634)
Reuses a pre-allocated matrix in tests to avoid allocating a fresh matrix in every call.
1 parent bb35dc9 commit 53048b2

File tree

1 file changed

+108
-107
lines changed

1 file changed

+108
-107
lines changed

stdlib/LinearAlgebra/test/triangular.jl

+108-107
Original file line numberDiff line numberDiff line change
@@ -35,6 +35,7 @@ debug && println("Test basic type functionality")
3535

3636
# Construct test matrix
3737
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')))
38+
M1 = Matrix(A1)
3839
@test t1(A1) === A1
3940
@test t1{elty1}(A1) === A1
4041
# test the ctor works for AbstractMatrix
@@ -68,7 +69,7 @@ debug && println("Test basic type functionality")
6869
@test simA1 == A1
6970

7071
# getindex
71-
let mA1 = Matrix(A1)
72+
let mA1 = M1
7273
# linear indexing
7374
for i in 1:length(A1)
7475
@test A1[i] == mA1[i]
@@ -141,8 +142,8 @@ debug && println("Test basic type functionality")
141142
#tril/triu
142143
if uplo1 === :L
143144
@test tril(A1,0) == A1
144-
@test tril(A1,-1) == LowerTriangular(tril(Matrix(A1), -1))
145-
@test tril(A1,1) == t1(tril(tril(Matrix(A1), 1)))
145+
@test tril(A1,-1) == LowerTriangular(tril(M1, -1))
146+
@test tril(A1,1) == t1(tril(tril(M1, 1)))
146147
@test tril(A1, -n - 2) == zeros(size(A1))
147148
@test tril(A1, n) == A1
148149
@test triu(A1,0) == t1(diagm(0 => diag(A1)))
@@ -152,8 +153,8 @@ debug && println("Test basic type functionality")
152153
@test triu(A1, n + 2) == zeros(size(A1))
153154
else
154155
@test triu(A1,0) == A1
155-
@test triu(A1,1) == UpperTriangular(triu(Matrix(A1), 1))
156-
@test triu(A1,-1) == t1(triu(triu(Matrix(A1), -1)))
156+
@test triu(A1,1) == UpperTriangular(triu(M1, 1))
157+
@test triu(A1,-1) == t1(triu(triu(M1, -1)))
157158
@test triu(A1, -n) == A1
158159
@test triu(A1, n + 2) == zeros(size(A1))
159160
@test tril(A1,0) == t1(diagm(0 => diag(A1)))
@@ -169,10 +170,10 @@ debug && println("Test basic type functionality")
169170
# [c]transpose[!] (test views as well, see issue #14317)
170171
let vrange = 1:n-1, viewA1 = t1(view(A1.data, vrange, vrange))
171172
# transpose
172-
@test copy(transpose(A1)) == transpose(Matrix(A1))
173+
@test copy(transpose(A1)) == transpose(M1)
173174
@test copy(transpose(viewA1)) == transpose(Matrix(viewA1))
174175
# adjoint
175-
@test copy(A1') == Matrix(A1)'
176+
@test copy(A1') == M1'
176177
@test copy(viewA1') == Matrix(viewA1)'
177178
# transpose!
178179
@test transpose!(copy(A1)) == transpose(A1)
@@ -185,28 +186,28 @@ debug && println("Test basic type functionality")
185186
end
186187

187188
# diag
188-
@test diag(A1) == diag(Matrix(A1))
189+
@test diag(A1) == diag(M1)
189190

190191
# tr
191-
@test tr(A1)::elty1 == tr(Matrix(A1))
192+
@test tr(A1)::elty1 == tr(M1)
192193

193194
# real
194-
@test real(A1) == real(Matrix(A1))
195-
@test imag(A1) == imag(Matrix(A1))
196-
@test abs.(A1) == abs.(Matrix(A1))
195+
@test real(A1) == real(M1)
196+
@test imag(A1) == imag(M1)
197+
@test abs.(A1) == abs.(M1)
197198

198199
# zero
199200
if A1 isa UpperTriangular || A1 isa LowerTriangular
200201
@test zero(A1) == zero(parent(A1))
201202
end
202203

203204
# Unary operations
204-
@test -A1 == -Matrix(A1)
205+
@test -A1 == -M1
205206

206207
# copy and copyto! (test views as well, see issue #14317)
207208
let vrange = 1:n-1, viewA1 = t1(view(A1.data, vrange, vrange))
208209
# copy
209-
@test copy(A1) == copy(Matrix(A1))
210+
@test copy(A1) == copy(M1)
210211
@test copy(viewA1) == copy(Matrix(viewA1))
211212
# copyto!
212213
B = similar(A1)
@@ -283,25 +284,25 @@ debug && println("Test basic type functionality")
283284
end
284285

285286
# Binary operations
286-
@test A1*0.5 == Matrix(A1)*0.5
287-
@test 0.5*A1 == 0.5*Matrix(A1)
288-
@test A1/0.5 == Matrix(A1)/0.5
289-
@test 0.5\A1 == 0.5\Matrix(A1)
287+
@test A1*0.5 == M1*0.5
288+
@test 0.5*A1 == 0.5*M1
289+
@test A1/0.5 == M1/0.5
290+
@test 0.5\A1 == 0.5\M1
290291

291292
# inversion
292-
@test inv(A1) inv(lu(Matrix(A1)))
293-
inv(Matrix(A1)) # issue #11298
293+
@test inv(A1) inv(lu(M1))
294+
inv(M1) # issue #11298
294295
@test isa(inv(A1), t1)
295296
# make sure the call to LAPACK works right
296297
if elty1 <: BlasFloat
297-
@test LinearAlgebra.inv!(copy(A1)) inv(lu(Matrix(A1)))
298+
@test LinearAlgebra.inv!(copy(A1)) inv(lu(M1))
298299
end
299300

300301
# Determinant
301-
@test det(A1) det(lu(Matrix(A1))) atol=sqrt(eps(real(float(one(elty1)))))*n*n
302-
@test logdet(A1) logdet(lu(Matrix(A1))) atol=sqrt(eps(real(float(one(elty1)))))*n*n
302+
@test det(A1) det(lu(M1)) atol=sqrt(eps(real(float(one(elty1)))))*n*n
303+
@test logdet(A1) logdet(lu(M1)) atol=sqrt(eps(real(float(one(elty1)))))*n*n
303304
lada, ladb = logabsdet(A1)
304-
flada, fladb = logabsdet(lu(Matrix(A1)))
305+
flada, fladb = logabsdet(lu(M1))
305306
@test lada flada atol=sqrt(eps(real(float(one(elty1)))))*n*n
306307
@test ladb fladb atol=sqrt(eps(real(float(one(elty1)))))*n*n
307308

@@ -324,7 +325,7 @@ debug && println("Test basic type functionality")
324325
for p in (1.0, Inf)
325326
@test cond(A1,p) cond(A1,p) atol=(cond(A1,p)+cond(A1,p))
326327
end
327-
@test cond(A1,2) == cond(Matrix(A1),2)
328+
@test cond(A1,2) == cond(M1,2)
328329
end
329330

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

336-
@test ((A1*A1)::t1) Matrix(A1) * Matrix(A1)
337-
@test ((A1/A1)::t1) Matrix(A1) / Matrix(A1)
338-
@test ((A1\A1)::t1) Matrix(A1) \ Matrix(A1)
337+
@test ((A1*A1)::t1) M1 * M1
338+
@test ((A1/A1)::t1) M1 / M1
339+
@test ((A1\A1)::t1) M1 \ M1
339340

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

349350
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')))
350-
351+
M2 = Matrix(A2)
351352
# Convert
352353
if elty1 <: Real && !(elty2 <: Integer)
353354
@test convert(AbstractMatrix{elty2}, A1) == t1(convert(Matrix{elty2}, A1.data))
@@ -356,21 +357,21 @@ debug && println("Test basic type functionality")
356357
end
357358

358359
# Binary operations
359-
@test A1 + A2 == Matrix(A1) + Matrix(A2)
360-
@test A1 - A2 == Matrix(A1) - Matrix(A2)
360+
@test A1 + A2 == M1 + M2
361+
@test A1 - A2 == M1 - M2
361362

362363
# Triangular-Triangular multiplication and division
363-
@test A1*A2 Matrix(A1)*Matrix(A2)
364-
@test transpose(A1)*A2 transpose(Matrix(A1))*Matrix(A2)
365-
@test transpose(A1)*adjoint(A2) transpose(Matrix(A1))*adjoint(Matrix(A2))
366-
@test adjoint(A1)*transpose(A2) adjoint(Matrix(A1))*transpose(Matrix(A2))
367-
@test A1'A2 Matrix(A1)'Matrix(A2)
368-
@test A1*transpose(A2) Matrix(A1)*transpose(Matrix(A2))
369-
@test A1*A2' Matrix(A1)*Matrix(A2)'
370-
@test transpose(A1)*transpose(A2) transpose(Matrix(A1))*transpose(Matrix(A2))
371-
@test A1'A2' Matrix(A1)'Matrix(A2)'
372-
@test A1/A2 Matrix(A1)/Matrix(A2)
373-
@test A1\A2 Matrix(A1)\Matrix(A2)
364+
@test A1*A2 M1*M2
365+
@test transpose(A1)*A2 transpose(M1)*M2
366+
@test transpose(A1)*adjoint(A2) transpose(M1)*adjoint(M2)
367+
@test adjoint(A1)*transpose(A2) adjoint(M1)*transpose(M2)
368+
@test A1'A2 M1'M2
369+
@test A1*transpose(A2) M1*transpose(M2)
370+
@test A1*A2' M1*M2'
371+
@test transpose(A1)*transpose(A2) transpose(M1)*transpose(M2)
372+
@test A1'A2' M1'M2'
373+
@test A1/A2 M1/M2
374+
@test A1\A2 M1\M2
374375
if uplo1 === :U && uplo2 === :U
375376
if t1 === UnitUpperTriangular && t2 === UnitUpperTriangular
376377
@test A1*A2 isa UnitUpperTriangular
@@ -411,20 +412,20 @@ debug && println("Test basic type functionality")
411412
@test_throws DimensionMismatch A2' * offsizeA
412413
@test_throws DimensionMismatch A2 * offsizeA
413414
if (uplo1 == uplo2 && elty1 == elty2 != Int && t1 != UnitLowerTriangular && t1 != UnitUpperTriangular)
414-
@test rdiv!(copy(A1), A2)::t1 A1/A2 Matrix(A1)/Matrix(A2)
415-
@test ldiv!(A2, copy(A1))::t1 A2\A1 Matrix(A2)\Matrix(A1)
415+
@test rdiv!(copy(A1), A2)::t1 A1/A2 M1/M2
416+
@test ldiv!(A2, copy(A1))::t1 A2\A1 M2\M1
416417
end
417418
if (uplo1 != uplo2 && elty1 == elty2 != Int && t2 != UnitLowerTriangular && t2 != UnitUpperTriangular)
418-
@test lmul!(adjoint(A1), copy(A2)) A1'*A2 Matrix(A1)'*Matrix(A2)
419-
@test lmul!(transpose(A1), copy(A2)) transpose(A1)*A2 transpose(Matrix(A1))*Matrix(A2)
420-
@test ldiv!(adjoint(A1), copy(A2)) A1'\A2 Matrix(A1)'\Matrix(A2)
421-
@test ldiv!(transpose(A1), copy(A2)) transpose(A1)\A2 transpose(Matrix(A1))\Matrix(A2)
419+
@test lmul!(adjoint(A1), copy(A2)) A1'*A2 M1'*M2
420+
@test lmul!(transpose(A1), copy(A2)) transpose(A1)*A2 transpose(M1)*M2
421+
@test ldiv!(adjoint(A1), copy(A2)) A1'\A2 M1'\M2
422+
@test ldiv!(transpose(A1), copy(A2)) transpose(A1)\A2 transpose(M1)\M2
422423
end
423424
if (uplo1 != uplo2 && elty1 == elty2 != Int && t1 != UnitLowerTriangular && t1 != UnitUpperTriangular)
424-
@test rmul!(copy(A1), adjoint(A2)) A1*A2' Matrix(A1)*Matrix(A2)'
425-
@test rmul!(copy(A1), transpose(A2)) A1*transpose(A2) Matrix(A1)*transpose(Matrix(A2))
426-
@test rdiv!(copy(A1), adjoint(A2)) A1/A2' Matrix(A1)/Matrix(A2)'
427-
@test rdiv!(copy(A1), transpose(A2)) A1/transpose(A2) Matrix(A1)/transpose(Matrix(A2))
425+
@test rmul!(copy(A1), adjoint(A2)) A1*A2' M1*M2'
426+
@test rmul!(copy(A1), transpose(A2)) A1*transpose(A2) M1*transpose(M2)
427+
@test rdiv!(copy(A1), adjoint(A2)) A1/A2' M1/M2'
428+
@test rdiv!(copy(A1), transpose(A2)) A1/transpose(A2) M1/transpose(M2)
428429
end
429430
end
430431
end
@@ -435,55 +436,55 @@ debug && println("Test basic type functionality")
435436
debug && println("elty1: $elty1, A1: $t1, B: $eltyB")
436437

437438
Tri = Tridiagonal(rand(eltyB,n-1),rand(eltyB,n),rand(eltyB,n-1))
438-
@test lmul!(Tri,copy(A1)) Tri*Matrix(A1)
439+
@test lmul!(Tri,copy(A1)) Tri*M1
439440
Tri = Tridiagonal(rand(eltyB,n-1),rand(eltyB,n),rand(eltyB,n-1))
440441
C = Matrix{promote_type(elty1,eltyB)}(undef, n, n)
441442
mul!(C, Tri, A1)
442-
@test C Tri*Matrix(A1)
443+
@test C Tri*M1
443444
Tri = Tridiagonal(rand(eltyB,n-1),rand(eltyB,n),rand(eltyB,n-1))
444445
mul!(C, A1, Tri)
445-
@test C Matrix(A1)*Tri
446+
@test C M1*Tri
446447

447448
# Triangular-dense Matrix/vector multiplication
448-
@test A1*B[:,1] Matrix(A1)*B[:,1]
449-
@test A1*B Matrix(A1)*B
450-
@test transpose(A1)*B[:,1] transpose(Matrix(A1))*B[:,1]
451-
@test A1'B[:,1] Matrix(A1)'B[:,1]
452-
@test transpose(A1)*B transpose(Matrix(A1))*B
453-
@test A1'B Matrix(A1)'B
454-
@test A1*transpose(B) Matrix(A1)*transpose(B)
455-
@test adjoint(A1)*transpose(B) Matrix(A1)'*transpose(B)
456-
@test transpose(A1)*adjoint(B) transpose(Matrix(A1))*adjoint(B)
457-
@test A1*B' Matrix(A1)*B'
458-
@test B*A1 B*Matrix(A1)
459-
@test transpose(B[:,1])*A1 transpose(B[:,1])*Matrix(A1)
460-
@test B[:,1]'A1 B[:,1]'Matrix(A1)
461-
@test transpose(B)*A1 transpose(B)*Matrix(A1)
462-
@test transpose(B)*adjoint(A1) transpose(B)*Matrix(A1)'
463-
@test adjoint(B)*transpose(A1) adjoint(B)*transpose(Matrix(A1))
464-
@test B'A1 B'Matrix(A1)
465-
@test B*transpose(A1) B*transpose(Matrix(A1))
466-
@test B*A1' B*Matrix(A1)'
467-
@test transpose(B[:,1])*transpose(A1) transpose(B[:,1])*transpose(Matrix(A1))
468-
@test B[:,1]'A1' B[:,1]'Matrix(A1)'
469-
@test transpose(B)*transpose(A1) transpose(B)*transpose(Matrix(A1))
470-
@test B'A1' B'Matrix(A1)'
449+
@test A1*B[:,1] M1*B[:,1]
450+
@test A1*B M1*B
451+
@test transpose(A1)*B[:,1] transpose(M1)*B[:,1]
452+
@test A1'B[:,1] M1'B[:,1]
453+
@test transpose(A1)*B transpose(M1)*B
454+
@test A1'B M1'B
455+
@test A1*transpose(B) M1*transpose(B)
456+
@test adjoint(A1)*transpose(B) M1'*transpose(B)
457+
@test transpose(A1)*adjoint(B) transpose(M1)*adjoint(B)
458+
@test A1*B' M1*B'
459+
@test B*A1 B*M1
460+
@test transpose(B[:,1])*A1 transpose(B[:,1])*M1
461+
@test B[:,1]'A1 B[:,1]'M1
462+
@test transpose(B)*A1 transpose(B)*M1
463+
@test transpose(B)*adjoint(A1) transpose(B)*M1'
464+
@test adjoint(B)*transpose(A1) adjoint(B)*transpose(M1)
465+
@test B'A1 B'M1
466+
@test B*transpose(A1) B*transpose(M1)
467+
@test B*A1' B*M1'
468+
@test transpose(B[:,1])*transpose(A1) transpose(B[:,1])*transpose(M1)
469+
@test B[:,1]'A1' B[:,1]'M1'
470+
@test transpose(B)*transpose(A1) transpose(B)*transpose(M1)
471+
@test B'A1' B'M1'
471472

472473
if eltyB == elty1
473-
@test mul!(similar(B), A1, B) Matrix(A1)*B
474-
@test mul!(similar(B), A1, adjoint(B)) Matrix(A1)*B'
475-
@test mul!(similar(B), A1, transpose(B)) Matrix(A1)*transpose(B)
476-
@test mul!(similar(B), adjoint(A1), adjoint(B)) Matrix(A1)'*B'
477-
@test mul!(similar(B), transpose(A1), transpose(B)) transpose(Matrix(A1))*transpose(B)
478-
@test mul!(similar(B), transpose(A1), adjoint(B)) transpose(Matrix(A1))*B'
479-
@test mul!(similar(B), adjoint(A1), transpose(B)) Matrix(A1)'*transpose(B)
480-
@test mul!(similar(B), adjoint(A1), B) Matrix(A1)'*B
481-
@test mul!(similar(B), transpose(A1), B) transpose(Matrix(A1))*B
474+
@test mul!(similar(B), A1, B) M1*B
475+
@test mul!(similar(B), A1, adjoint(B)) M1*B'
476+
@test mul!(similar(B), A1, transpose(B)) M1*transpose(B)
477+
@test mul!(similar(B), adjoint(A1), adjoint(B)) M1'*B'
478+
@test mul!(similar(B), transpose(A1), transpose(B)) transpose(M1)*transpose(B)
479+
@test mul!(similar(B), transpose(A1), adjoint(B)) transpose(M1)*B'
480+
@test mul!(similar(B), adjoint(A1), transpose(B)) M1'*transpose(B)
481+
@test mul!(similar(B), adjoint(A1), B) M1'*B
482+
@test mul!(similar(B), transpose(A1), B) transpose(M1)*B
482483
# test also vector methods
483484
B1 = vec(B[1,:])
484-
@test mul!(similar(B1), A1, B1) Matrix(A1)*B1
485-
@test mul!(similar(B1), adjoint(A1), B1) Matrix(A1)'*B1
486-
@test mul!(similar(B1), transpose(A1), B1) transpose(Matrix(A1))*B1
485+
@test mul!(similar(B1), A1, B1) M1*B1
486+
@test mul!(similar(B1), adjoint(A1), B1) M1'*B1
487+
@test mul!(similar(B1), transpose(A1), B1) transpose(M1)*B1
487488
end
488489
#error handling
489490
Ann, Bmm, bm = A1, Matrix{eltyB}(undef, n+1, n+1), Vector{eltyB}(undef, n+1)
@@ -495,30 +496,30 @@ debug && println("Test basic type functionality")
495496
@test_throws DimensionMismatch rmul!(Bmm, transpose(Ann))
496497

497498
# ... and division
498-
@test A1\B[:,1] Matrix(A1)\B[:,1]
499-
@test A1\B Matrix(A1)\B
500-
@test transpose(A1)\B[:,1] transpose(Matrix(A1))\B[:,1]
501-
@test A1'\B[:,1] Matrix(A1)'\B[:,1]
502-
@test transpose(A1)\B transpose(Matrix(A1))\B
503-
@test A1'\B Matrix(A1)'\B
504-
@test A1\transpose(B) Matrix(A1)\transpose(B)
505-
@test A1\B' Matrix(A1)\B'
506-
@test transpose(A1)\transpose(B) transpose(Matrix(A1))\transpose(B)
507-
@test A1'\B' Matrix(A1)'\B'
499+
@test A1\B[:,1] M1\B[:,1]
500+
@test A1\B M1\B
501+
@test transpose(A1)\B[:,1] transpose(M1)\B[:,1]
502+
@test A1'\B[:,1] M1'\B[:,1]
503+
@test transpose(A1)\B transpose(M1)\B
504+
@test A1'\B M1'\B
505+
@test A1\transpose(B) M1\transpose(B)
506+
@test A1\B' M1\B'
507+
@test transpose(A1)\transpose(B) transpose(M1)\transpose(B)
508+
@test A1'\B' M1'\B'
508509
Ann, bm = A1, Vector{elty1}(undef,n+1)
509510
@test_throws DimensionMismatch Ann\bm
510511
@test_throws DimensionMismatch Ann'\bm
511512
@test_throws DimensionMismatch transpose(Ann)\bm
512513
if t1 == UpperTriangular || t1 == LowerTriangular
513514
@test_throws SingularException ldiv!(t1(zeros(elty1, n, n)), fill(eltyB(1), n))
514515
end
515-
@test B/A1 B/Matrix(A1)
516-
@test B/transpose(A1) B/transpose(Matrix(A1))
517-
@test B/A1' B/Matrix(A1)'
518-
@test transpose(B)/A1 transpose(B)/Matrix(A1)
519-
@test B'/A1 B'/Matrix(A1)
520-
@test transpose(B)/transpose(A1) transpose(B)/transpose(Matrix(A1))
521-
@test B'/A1' B'/Matrix(A1)'
516+
@test B/A1 B/M1
517+
@test B/transpose(A1) B/transpose(M1)
518+
@test B/A1' B/M1'
519+
@test transpose(B)/A1 transpose(B)/M1
520+
@test B'/A1 B'/M1
521+
@test transpose(B)/transpose(A1) transpose(B)/transpose(M1)
522+
@test B'/A1' B'/M1'
522523

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

0 commit comments

Comments
 (0)