From 80ec7925cf1803196750916b1f7ca6a4700c112d Mon Sep 17 00:00:00 2001 From: Alexis Montoison Date: Tue, 14 May 2024 12:26:05 +0200 Subject: [PATCH] Update test_BLAS.jl --- test/test_BLAS.jl | 139 ++++++++++++++++++++++++++++------------------ 1 file changed, 84 insertions(+), 55 deletions(-) diff --git a/test/test_BLAS.jl b/test/test_BLAS.jl index c27db96..f6eb864 100644 --- a/test/test_BLAS.jl +++ b/test/test_BLAS.jl @@ -25,74 +25,103 @@ macro test_blas(ex) end @testset "matrix-vector multiplication (non-square)" begin - for i = 1:5 - a = sprand(10, 5, 0.5) - b = rand(5) - @test_blas maximum(abs.(a*b - Array(a)*b)) < 100*eps() - b = rand(5, 5) - @test_blas maximum(abs.(a*b - Array(a)*b)) < 100*eps() - b = rand(10) - @test_blas maximum(abs.(a'*b - Array(a)'*b)) < 100*eps() - @test_blas maximum(abs.(transpose(a)*b - Array(a)'*b)) < 100*eps() - b = rand(10,10) - @test_blas maximum(abs.(a'*b - Array(a)'*b)) < 100*eps() - @test_blas maximum(abs.(transpose(a)*b - Array(a)'*b)) < 100*eps() + T = Float64 + @testset "interface = $interface" for interface in ("LP64", "ILP64") + S = interface == "LP64" ? Int32 : Int64 + for i = 1:5 + a = sprand(T, 10, 5, 0.5) + a = SparseMatrixCSC{T, S}(a) + b = rand(T, 5) + @test_blas maximum(abs.(a*b - Array(a)*b)) < 100*eps(T) + b = rand(T, 5, 5) + @test_blas maximum(abs.(a*b - Array(a)*b)) < 100*eps(T) + b = rand(T, 10) + @test_blas maximum(abs.(a'*b - Array(a)'*b)) < 100*eps(T) + @test_blas maximum(abs.(transpose(a)*b - Array(a)'*b)) < 100*eps(T) + b = rand(T, 10, 10) + @test_blas maximum(abs.(a'*b - Array(a)'*b)) < 100*eps() + @test_blas maximum(abs.(transpose(a)*b - Array(a)'*b)) < 100*eps(T) + end end end #? @testset "complex matrix-vector multiplication" begin - for i = 1:5 - a = I + im * 0.1*sprandn(5, 5, 0.2) - b = randn(5,3) + im*randn(5,3) - c = randn(5) + im*randn(5) - d = randn(5) + im*randn(5) - α = rand(ComplexF64) - β = rand(ComplexF64) - @test_blas (maximum(abs.(a*b - Array(a)*b)) < 100*eps()) - @test_blas (maximum(abs.(a'*b - Array(a)'*b)) < 100*eps()) - @test_blas (maximum(abs.(transpose(a)*b - transpose(Array(a))*b)) < 100*eps()) - @test_blas (maximum(abs.(mul!(similar(b), a, b) - Array(a)*b)) < 100*eps()) - @test_blas (maximum(abs.(mul!(similar(c), a, c) - Array(a)*c)) < 100*eps()) - @test_blas (maximum(abs.(mul!(similar(b), transpose(a), b) - transpose(Array(a))*b)) < 100*eps()) - @test_blas (maximum(abs.(mul!(similar(c), transpose(a), c) - transpose(Array(a))*c)) < 100*eps()) - @test_blas (maximum(abs.(mul!(copy(b), a, b, α, β) - (α*(Array(a)*b) + β*b))) < 100*eps()) - @test_blas (maximum(abs.(mul!(copy(b), transpose(a), b, α, β) - (α*(transpose(Array(a))*b) + β*b))) < 100*eps()) - @test_blas (maximum(abs.(mul!(copy(c), transpose(a), c, α, β) - (α*(transpose(Array(a))*c) + β*c))) < 100*eps()) - α = β = 1 # test conversion to float - @test_blas (maximum(abs.(mul!(copy(b), a, b, α, β) - (α*(Array(a)*b) + β*b))) < 100*eps()) - @test_blas (maximum(abs.(mul!(copy(b), transpose(a), b, α, β) - (α*(transpose(Array(a))*b) + β*b))) < 100*eps()) - @test_blas (maximum(abs.(mul!(copy(c), transpose(a), c, α, β) - (α*(transpose(Array(a))*c) + β*c))) < 100*eps()) + T = ComplexF64 + R = Float64 + @testset "interface = $interface" for interface in ("LP64", "ILP64") + S = interface == "LP64" ? Int32 : Int64 + for i = 1:5 + a = I + im * 0.1 * sprandn(T, 5, 5, 0.2) + a = SparseMatrixCSC{T,S}(a) + b = randn(T, 5, 3) + im * randn(T, 5, 3) + c = randn(T, 5) + im * randn(T, 5) + d = randn(T, 5) + im * randn(T, 5) + α = rand(T) + β = rand(T) + @test_blas (maximum(abs.(a*b - Array(a)*b)) < 100*eps(R)) + @test_blas (maximum(abs.(a'*b - Array(a)'*b)) < 100*eps(R)) + @test_blas (maximum(abs.(transpose(a)*b - transpose(Array(a))*b)) < 100*eps(R)) + @test_blas (maximum(abs.(mul!(similar(b), a, b) - Array(a)*b)) < 100*eps(R)) + @test_blas (maximum(abs.(mul!(similar(c), a, c) - Array(a)*c)) < 100*eps(R)) + @test_blas (maximum(abs.(mul!(similar(b), transpose(a), b) - transpose(Array(a))*b)) < 100*eps(R)) + @test_blas (maximum(abs.(mul!(similar(c), transpose(a), c) - transpose(Array(a))*c)) < 100*eps(R)) + @test_blas (maximum(abs.(mul!(copy(b), a, b, α, β) - (α*(Array(a)*b) + β*b))) < 100*eps(R)) + @test_blas (maximum(abs.(mul!(copy(b), transpose(a), b, α, β) - (α*(transpose(Array(a))*b) + β*b))) < 100*eps(R)) + @test_blas (maximum(abs.(mul!(copy(c), transpose(a), c, α, β) - (α*(transpose(Array(a))*c) + β*c))) < 100*eps(R)) + α = β = 1 # test conversion to float + @test_blas (maximum(abs.(mul!(copy(b), a, b, α, β) - (α*(Array(a)*b) + β*b))) < 100*eps(R)) + @test_blas (maximum(abs.(mul!(copy(b), transpose(a), b, α, β) - (α*(transpose(Array(a))*b) + β*b))) < 100*eps(R)) + @test_blas (maximum(abs.(mul!(copy(c), transpose(a), c, α, β) - (α*(transpose(Array(a))*c) + β*c))) < 100*eps(R)) - c = randn(6) + im*randn(6) - @test_throws DimensionMismatch transpose(a)*c - @test_throws DimensionMismatch a.*c - @test_throws DimensionMismatch a.*c + c = randn(T, 6) + im * randn(T, 6) + @test_throws DimensionMismatch transpose(a)*c + @test_throws DimensionMismatch a.*c + @test_throws DimensionMismatch a.*c + end end end @testset "triangular" begin - n = 100 - A = sprandn(n, n, 0.5) + sqrt(n)*I - b = rand(n) - symA = A + transpose(A) - trilA = tril(A) - triuA = triu(A) - trilUA = tril(A, -1) + I - triuUA = triu(A, 1) + I + @testset "T = $T" for T in (Float32, Float64, ComplexF32, ComplexF64) + @testset "interface = $interface" for interface in ("LP64", "ILP64") + S = interface == "LP64" ? Int32 : Int64 + n = 100 + A = sprandn(T, n, n, 0.5) + sqrt(n)*I + A = SparseMatrixCSC{T,S}(A) + b = rand(T, n) + trilA = tril(A) + triuA = triu(A) + trilUA = tril(A, -1) + I + triuUA = triu(A, 1) + I - @test_blas LowerTriangular(trilA) \ b ≈ Array(LowerTriangular(trilA)) \ b - @test_blas LowerTriangular(trilA) * b ≈ Array(LowerTriangular(trilA)) * b + @test_blas LowerTriangular(trilA) \ b ≈ Array(LowerTriangular(trilA)) \ b + @test_blas LowerTriangular(trilA) * b ≈ Array(LowerTriangular(trilA)) * b - @test_blas UpperTriangular(triuA) \ b ≈ Array(UpperTriangular(triuA)) \ b - @test_blas UpperTriangular(triuA) * b ≈ Array(UpperTriangular(triuA)) * b + @test_blas UpperTriangular(triuA) \ b ≈ Array(UpperTriangular(triuA)) \ b + @test_blas UpperTriangular(triuA) * b ≈ Array(UpperTriangular(triuA)) * b - @test_blas UnitLowerTriangular(trilUA) \ b ≈ Array(UnitLowerTriangular(trilUA)) \ b - @test_blas UnitLowerTriangular(trilUA) * b ≈ Array(UnitLowerTriangular(trilUA)) * b + @test_blas UnitLowerTriangular(trilUA) \ b ≈ Array(UnitLowerTriangular(trilUA)) \ b + @test_blas UnitLowerTriangular(trilUA) * b ≈ Array(UnitLowerTriangular(trilUA)) * b - @test_blas UnitUpperTriangular(triuUA) \ b ≈ Array(UnitUpperTriangular(triuUA)) \ b - @test_blas UnitUpperTriangular(triuUA) * b ≈ Array(UnitUpperTriangular(triuUA)) * b + @test_blas UnitUpperTriangular(triuUA) \ b ≈ Array(UnitUpperTriangular(triuUA)) \ b + @test_blas UnitUpperTriangular(triuUA) * b ≈ Array(UnitUpperTriangular(triuUA)) * b + end + end +end - @test_blas Symmetric(symA) * b ≈ Array(Symmetric(symA)) * b - @test_blas Hermitian(symA) * b ≈ Array(Hermitian(symA)) * b +@testset "Symmetric -- Hermitian" begin + @testset "T = $T" for T in (Float32, Float64, ComplexF32, ComplexF64) + @testset "interface = $interface" for interface in ("LP64", "ILP64") + S = interface == "LP64" ? Int32 : Int64 + n = 100 + A = sprandn(T, n, n, 0.5) + sqrt(n)*I + b = rand(T, n) + symA = A + transpose(A) + hermA = A + adjoint(A) + + @test_blas Symmetric(symA) * b ≈ Array(Symmetric(symA)) * b + @test_blas Hermitian(hermA) * b ≈ Array(Hermitian(hermA)) * b + end + end end