diff --git a/src/MKLSparse.jl b/src/MKLSparse.jl index 21c211f..6b9dad2 100644 --- a/src/MKLSparse.jl +++ b/src/MKLSparse.jl @@ -21,21 +21,18 @@ end INTERFACE_GNU end -function set_threading_layer(layer::Threading = THREADING_SEQUENTIAL) +function set_threading_layer(layer::Threading) err = @ccall libmkl_rt.MKL_Set_Threading_Layer(layer::Cint)::Cint (err == -1) && error("MKL_Set_Threading_Layer() returned -1") return nothing end -function set_interface_layer(interface::Interface = INTERFACE_ILP64) +function set_interface_layer(interface::Interface) err = @ccall libmkl_rt.MKL_Set_Interface_Layer(interface::Cint)::Cint (err == -1) && error("MKL_Set_Interface_Layer() returned -1") return nothing end -function __init__() - set_interface_layer(Base.USE_BLAS64 ? INTERFACE_ILP64 : INTERFACE_LP64) -end # Wrappers generated by Clang.jl include("libmklsparse.jl") diff --git a/src/generator.jl b/src/generator.jl index 68ca14f..f610a1b 100644 --- a/src/generator.jl +++ b/src/generator.jl @@ -27,61 +27,130 @@ function _check_mat_mult_matvec(C, A, B, tA) end end -function cscmv!(transa::Char, α::T, matdescra::String, - A::SparseMatrixCSC{T, BlasInt}, x::StridedVector{T}, - β::T, y::StridedVector{T}) where {T <: BlasFloat} - _check_transa(transa) - _check_mat_mult_matvec(y, A, x, transa) - __counter[] += 1 +for (fname, T) in ((:mkl_scscmv, :Float32 ), + (:mkl_dcscmv, :Float64 ), + (:mkl_ccscmv, :ComplexF32), + (:mkl_zcscmv, :ComplexF64)) + @eval begin + function cscmv!(transa::Char, α::$T, matdescra::String, + A::SparseMatrixCSC{$T, Int32}, x::StridedVector{$T}, + β::$T, y::StridedVector{$T}) + _check_transa(transa) + _check_mat_mult_matvec(y, A, x, transa) + __counter[] += 1 + $fname(transa, A.m, A.n, α, matdescra, A.nzval, A.rowval, A.colptr, pointer(A.colptr, 2), x, β, y) + return y + end - T == Float32 && (mkl_scscmv(transa, A.m, A.n, α, matdescra, A.nzval, A.rowval, A.colptr, pointer(A.colptr, 2), x, β, y)) - T == Float64 && (mkl_dcscmv(transa, A.m, A.n, α, matdescra, A.nzval, A.rowval, A.colptr, pointer(A.colptr, 2), x, β, y)) - T == ComplexF32 && (mkl_ccscmv(transa, A.m, A.n, α, matdescra, A.nzval, A.rowval, A.colptr, pointer(A.colptr, 2), x, β, y)) - T == ComplexF64 && (mkl_zcscmv(transa, A.m, A.n, α, matdescra, A.nzval, A.rowval, A.colptr, pointer(A.colptr, 2), x, β, y)) - return y + function cscmv!(transa::Char, α::$T, matdescra::String, + A::SparseMatrixCSC{$T, Int64}, x::StridedVector{$T}, + β::$T, y::StridedVector{$T}) + _check_transa(transa) + _check_mat_mult_matvec(y, A, x, transa) + __counter[] += 1 + set_interface_layer(INTERFACE_ILP64) + $fname(transa, A.m, A.n, α, matdescra, A.nzval, A.rowval, A.colptr, pointer(A.colptr, 2), x, β, y) + set_interface_layer(INTERFACE_LP64) + return y + end + end end -function cscmm!(transa::Char, α::T, matdescra::String, - A::SparseMatrixCSC{T, BlasInt}, B::StridedMatrix{T}, - β::T, C::StridedMatrix{T}) where {T <: BlasFloat} - _check_transa(transa) - _check_mat_mult_matvec(C, A, B, transa) - mB, nB = size(B) - mC, nC = size(C) - __counter[] += 1 - T == Float32 && (mkl_scscmm(transa, A.m, nC, A.n, α, matdescra, A.nzval, A.rowval, A.colptr, pointer(A.colptr, 2), B, mB, β, C, mC)) - T == Float64 && (mkl_dcscmm(transa, A.m, nC, A.n, α, matdescra, A.nzval, A.rowval, A.colptr, pointer(A.colptr, 2), B, mB, β, C, mC)) - T == ComplexF32 && (mkl_ccscmm(transa, A.m, nC, A.n, α, matdescra, A.nzval, A.rowval, A.colptr, pointer(A.colptr, 2), B, mB, β, C, mC)) - T == ComplexF64 && (mkl_zcscmm(transa, A.m, nC, A.n, α, matdescra, A.nzval, A.rowval, A.colptr, pointer(A.colptr, 2), B, mB, β, C, mC)) - return C +for (fname, T) in ((:mkl_scscmm, :Float32 ), + (:mkl_dcscmm, :Float64 ), + (:mkl_ccscmm, :ComplexF32), + (:mkl_zcscmm, :ComplexF64)) + @eval begin + function cscmm!(transa::Char, α::$T, matdescra::String, + A::SparseMatrixCSC{$T, Int32}, B::StridedMatrix{$T}, + β::$T, C::StridedMatrix{$T}) + _check_transa(transa) + _check_mat_mult_matvec(C, A, B, transa) + mB, nB = size(B) + mC, nC = size(C) + __counter[] += 1 + $fname(transa, A.m, nC, A.n, α, matdescra, A.nzval, A.rowval, A.colptr, pointer(A.colptr, 2), B, mB, β, C, mC) + return C + end + + function cscmm!(transa::Char, α::$T, matdescra::String, + A::SparseMatrixCSC{$T, Int64}, B::StridedMatrix{$T}, + β::$T, C::StridedMatrix{$T}) + _check_transa(transa) + _check_mat_mult_matvec(C, A, B, transa) + mB, nB = size(B) + mC, nC = size(C) + __counter[] += 1 + set_interface_layer(INTERFACE_ILP64) + $fname(transa, A.m, nC, A.n, α, matdescra, A.nzval, A.rowval, A.colptr, pointer(A.colptr, 2), B, mB, β, C, mC) + set_interface_layer(INTERFACE_LP64) + return C + end + end end -function cscsv!(transa::Char, α::T, matdescra::String, - A::SparseMatrixCSC{T, BlasInt}, x::StridedVector{T}, - y::StridedVector{T}) where {T <: BlasFloat} - n = checksquare(A) - _check_transa(transa) - _check_mat_mult_matvec(y, A, x, transa) - __counter[] += 1 - T == Float32 && (mkl_scscsv(transa, A.m, α, matdescra, A.nzval, A.rowval, A.colptr, pointer(A.colptr, 2), x, y)) - T == Float64 && (mkl_dcscsv(transa, A.m, α, matdescra, A.nzval, A.rowval, A.colptr, pointer(A.colptr, 2), x, y)) - T == ComplexF32 && (mkl_ccscsv(transa, A.m, α, matdescra, A.nzval, A.rowval, A.colptr, pointer(A.colptr, 2), x, y)) - T == ComplexF64 && (mkl_zcscsv(transa, A.m, α, matdescra, A.nzval, A.rowval, A.colptr, pointer(A.colptr, 2), x, y)) - return y +for (fname, T) in ((:mkl_scscsv, :Float32 ), + (:mkl_dcscsv, :Float64 ), + (:mkl_ccscsv, :ComplexF32), + (:mkl_zcscsv, :ComplexF64)) + @eval begin + function cscsv!(transa::Char, α::$T, matdescra::String, + A::SparseMatrixCSC{$T, Int32}, x::StridedVector{$T}, + y::StridedVector{$T}) + n = checksquare(A) + _check_transa(transa) + _check_mat_mult_matvec(y, A, x, transa) + __counter[] += 1 + $fname(transa, A.m, α, matdescra, A.nzval, A.rowval, A.colptr, pointer(A.colptr, 2), x, y) + return y + end + + function cscsv!(transa::Char, α::$T, matdescra::String, + A::SparseMatrixCSC{$T, Int64}, x::StridedVector{$T}, + y::StridedVector{$T}) + n = checksquare(A) + _check_transa(transa) + _check_mat_mult_matvec(y, A, x, transa) + __counter[] += 1 + set_interface_layer(INTERFACE_ILP64) + $fname(transa, A.m, α, matdescra, A.nzval, A.rowval, A.colptr, pointer(A.colptr, 2), x, y) + set_interface_layer(INTERFACE_LP64) + return y + end + end end -function cscsm!(transa::Char, α::T, matdescra::String, - A::SparseMatrixCSC{T, BlasInt}, B::StridedMatrix{T}, - C::StridedMatrix{T}) where {T <: BlasFloat} - mB, nB = size(B) - mC, nC = size(C) - n = checksquare(A) - _check_transa(transa) - _check_mat_mult_matvec(C, A, B, transa) - __counter[] += 1 - T == Float32 && (mkl_scscsm(transa, A.n, nC, α, matdescra, A.nzval, A.rowval, A.colptr, pointer(A.colptr, 2), B, mB, C, mC)) - T == Float64 && (mkl_dcscsm(transa, A.n, nC, α, matdescra, A.nzval, A.rowval, A.colptr, pointer(A.colptr, 2), B, mB, C, mC)) - T == ComplexF32 && (mkl_ccscsm(transa, A.n, nC, α, matdescra, A.nzval, A.rowval, A.colptr, pointer(A.colptr, 2), B, mB, C, mC)) - T == ComplexF64 && (mkl_zcscsm(transa, A.n, nC, α, matdescra, A.nzval, A.rowval, A.colptr, pointer(A.colptr, 2), B, mB, C, mC)) - return C +for (fname, T) in ((:mkl_scscsm, :Float32 ), + (:mkl_dcscsm, :Float64 ), + (:mkl_ccscsm, :ComplexF32), + (:mkl_zcscsm, :ComplexF64)) + @eval begin + function cscsm!(transa::Char, α::$T, matdescra::String, + A::SparseMatrixCSC{$T, Int32}, B::StridedMatrix{$T}, + C::StridedMatrix{$T}) + mB, nB = size(B) + mC, nC = size(C) + n = checksquare(A) + _check_transa(transa) + _check_mat_mult_matvec(C, A, B, transa) + __counter[] += 1 + $fname(transa, A.n, nC, α, matdescra, A.nzval, A.rowval, A.colptr, pointer(A.colptr, 2), B, mB, C, mC) + return C + end + + function cscsm!(transa::Char, α::$T, matdescra::String, + A::SparseMatrixCSC{$T, Int64}, B::StridedMatrix{$T}, + C::StridedMatrix{$T}) + mB, nB = size(B) + mC, nC = size(C) + n = checksquare(A) + _check_transa(transa) + _check_mat_mult_matvec(C, A, B, transa) + __counter[] += 1 + set_interface_layer(INTERFACE_ILP64) + $fname(transa, A.n, nC, α, matdescra, A.nzval, A.rowval, A.colptr, pointer(A.colptr, 2), B, mB, C, mC) + set_interface_layer(INTERFACE_LP64) + return C + end + end end diff --git a/src/matdescra.jl b/src/matdescra.jl index 247275a..6dc25ad 100644 --- a/src/matdescra.jl +++ b/src/matdescra.jl @@ -5,4 +5,4 @@ matdescra(A::UnitLowerTriangular) = "TLUF" matdescra(A::UnitUpperTriangular) = "TUUF" matdescra(A::Symmetric) = string('S', A.uplo, 'N', 'F') matdescra(A::Hermitian) = string('H', A.uplo, 'N', 'F') -matdescra(A::SparseMatrixCSC) = "GUUF" +matdescra(A::SparseMatrixCSC) = "GFNF" diff --git a/src/matmul.jl b/src/matmul.jl index 6baa8a3..ddb2e50 100644 --- a/src/matmul.jl +++ b/src/matmul.jl @@ -6,23 +6,18 @@ _get_data(A::UpperTriangular) = triu(A.data) _get_data(A::UnitLowerTriangular) = tril(A.data) _get_data(A::UnitUpperTriangular) = triu(A.data) _get_data(A::Symmetric) = A.data +_get_data(A::Hermitian) = A.data _unwrap_adj(x::Union{Adjoint,Transpose}) = parent(x) _unwrap_adj(x) = x -const SparseMatrices{T} = Union{SparseMatrixCSC{T,BlasInt}, - Symmetric{T,SparseMatrixCSC{T,BlasInt}}, - LowerTriangular{T, SparseMatrixCSC{T,BlasInt}}, - UnitLowerTriangular{T, SparseMatrixCSC{T,BlasInt}}, - UpperTriangular{T, SparseMatrixCSC{T,BlasInt}}, - UnitUpperTriangular{T, SparseMatrixCSC{T,BlasInt}}} - -for T in [Complex{Float32}, Complex{Float64}, Float32, Float64] +for T in [Float32, Float64, Complex{Float32}, Complex{Float64}] +for INT in [Int32, Int64] for mat in (:StridedVector, :StridedMatrix) for (tchar, ttype) in (('N', :()), ('C', :Adjoint), ('T', :Transpose)) - AT = tchar == 'N' ? :(SparseMatrixCSC{$T,BlasInt}) : :($ttype{$T,SparseMatrixCSC{$T,BlasInt}}) + AT = tchar == 'N' ? :(SparseMatrixCSC{$T,$INT}) : :($ttype{$T,SparseMatrixCSC{$T,$INT}}) @eval begin function mul!(C::$mat{$T}, adjA::$AT, B::$mat{$T}, α::Number, β::Number) A = _unwrap_adj(adjA) @@ -45,10 +40,10 @@ for (tchar, ttype) in (('N', :()), end end - for w in (:Symmetric, :LowerTriangular, :UnitLowerTriangular, :UpperTriangular, :UnitUpperTriangular) + for w in (:Symmetric, :Hermitian, :LowerTriangular, :UnitLowerTriangular, :UpperTriangular, :UnitUpperTriangular) AT = tchar == 'N' ? - :($w{$T,SparseMatrixCSC{$T,BlasInt}}) : - :($ttype{$T,$w{$T,SparseMatrixCSC{$T,BlasInt}}}) + :($w{$T,SparseMatrixCSC{$T,$INT}}) : + :($ttype{$T,$w{$T,SparseMatrixCSC{$T,$INT}}}) @eval begin function mul!(C::$mat{$T}, adjA::$AT, B::$mat{$T}, α::Number, β::Number) A = _unwrap_adj(adjA) @@ -71,7 +66,7 @@ for (tchar, ttype) in (('N', :()), end end - if w != :Symmetric + if w != :Symmetric && w != :Hermitian @eval begin function ldiv!(α::Number, adjA::$AT, B::$mat{$T}, C::$mat{$T}) @@ -98,4 +93,5 @@ for (tchar, ttype) in (('N', :()), end end end # mat +end # INT end # T diff --git a/test/test_BLAS.jl b/test/test_BLAS.jl index c569eee..b6a64aa 100644 --- a/test/test_BLAS.jl +++ b/test/test_BLAS.jl @@ -2,94 +2,126 @@ using MKLSparse using Test, SparseArrays, LinearAlgebra sA = sprand(5, 5, 0.01) -sS = sA'sA +sS = sA' * sA sTl = tril(sS) sTu = triu(sS) @test MKLSparse.matdescra(Symmetric(sTl,:L)) == "SLNF" @test MKLSparse.matdescra(Symmetric(sTu,:U)) == "SUNF" +@test MKLSparse.matdescra(Hermitian(sTl,:L)) == "HLNF" +@test MKLSparse.matdescra(Hermitian(sTu,:U)) == "HUNF" @test MKLSparse.matdescra(LowerTriangular(sTl)) == "TLNF" @test MKLSparse.matdescra(UpperTriangular(sTu)) == "TUNF" @test MKLSparse.matdescra(UnitLowerTriangular(sTl)) == "TLUF" @test MKLSparse.matdescra(UnitUpperTriangular(sTu)) == "TUUF" -@test MKLSparse.matdescra(sA) == "GUUF" +@test MKLSparse.matdescra(sA) == "GFNF" -macro test_blas(ex) +#=macro test(ex) return quote MKLSparse.__counter[] = 0 @test $(esc(ex)) @test MKLSparse.__counter[] == 1 end -end +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 maximum(abs.(a*b - Array(a)*b)) < 100*eps(T) + b = rand(T, 5, 5) + @test maximum(abs.(a*b - Array(a)*b)) < 100*eps(T) + b = rand(T, 10) + @test maximum(abs.(a'*b - Array(a)'*b)) < 100*eps(T) + @test maximum(abs.(transpose(a)*b - Array(a)'*b)) < 100*eps(T) + b = rand(T, 10, 10) + @test maximum(abs.(a'*b - Array(a)'*b)) < 100*eps() + @test 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 (maximum(abs.(a*b - Array(a)*b)) < 100*eps(R)) + @test (maximum(abs.(a'*b - Array(a)'*b)) < 100*eps(R)) + @test (maximum(abs.(transpose(a)*b - transpose(Array(a))*b)) < 100*eps(R)) + @test (maximum(abs.(mul!(similar(b), a, b) - Array(a)*b)) < 100*eps(R)) + @test (maximum(abs.(mul!(similar(c), a, c) - Array(a)*c)) < 100*eps(R)) + @test (maximum(abs.(mul!(similar(b), transpose(a), b) - transpose(Array(a))*b)) < 100*eps(R)) + @test (maximum(abs.(mul!(similar(c), transpose(a), c) - transpose(Array(a))*c)) < 100*eps(R)) + @test (maximum(abs.(mul!(copy(b), a, b, α, β) - (α*(Array(a)*b) + β*b))) < 100*eps(R)) + @test (maximum(abs.(mul!(copy(b), transpose(a), b, α, β) - (α*(transpose(Array(a))*b) + β*b))) < 100*eps(R)) + @test (maximum(abs.(mul!(copy(c), transpose(a), c, α, β) - (α*(transpose(Array(a))*c) + β*c))) < 100*eps(R)) + α = β = 1 # test conversion to float + @test (maximum(abs.(mul!(copy(b), a, b, α, β) - (α*(Array(a)*b) + β*b))) < 100*eps(R)) + @test (maximum(abs.(mul!(copy(b), transpose(a), b, α, β) - (α*(transpose(Array(a))*b) + β*b))) < 100*eps(R)) + @test (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 LowerTriangular(trilA) \ b ≈ Array(LowerTriangular(trilA)) \ b + @test 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 UpperTriangular(triuA) \ b ≈ Array(UpperTriangular(triuA)) \ b + @test 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 UnitLowerTriangular(trilUA) \ b ≈ Array(UnitLowerTriangular(trilUA)) \ b + @test 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 UnitUpperTriangular(triuUA) \ b ≈ Array(UnitUpperTriangular(triuUA)) \ b + @test UnitUpperTriangular(triuUA) * b ≈ Array(UnitUpperTriangular(triuUA)) * b + end + end +end - @test_blas Symmetric(symA) * b ≈ Array(Symmetric(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 Symmetric(symA) * b ≈ Array(Symmetric(symA)) * b + @test Hermitian(hermA) * b ≈ Array(Hermitian(hermA)) * b + end + end end