Skip to content

Commit

Permalink
Merge pull request #16 from JuliaSparse/an/07
Browse files Browse the repository at this point in the history
Updates for Julia 0.7
  • Loading branch information
andreasnoack authored Jan 10, 2019
2 parents ce41674 + 4d63e62 commit 3957a66
Show file tree
Hide file tree
Showing 6 changed files with 126 additions and 110 deletions.
3 changes: 1 addition & 2 deletions REQUIRE
Original file line number Diff line number Diff line change
@@ -1,2 +1 @@
julia 0.5
Compat 0.18
julia 0.7
5 changes: 2 additions & 3 deletions src/BLAS/BLAS.jl
Original file line number Diff line number Diff line change
@@ -1,8 +1,7 @@
module BLAS

using Compat

import Base.LinAlg: checksquare, UnitLowerTriangular, UnitUpperTriangular, BlasInt, BlasFloat
using LinearAlgebra, SparseArrays
using LinearAlgebra: BlasInt, checksquare

# For testing purposes:
global const __counter = Ref(0)
Expand Down
58 changes: 29 additions & 29 deletions src/BLAS/level_2_3/generator.jl
Original file line number Diff line number Diff line change
Expand Up @@ -29,22 +29,22 @@ end

for (mv, sv, mm, sm, T) in ((:mkl_scscmv, :mkl_scscsv, :mkl_scscmm, :mkl_scscsm, :Float32),
(:mkl_dcscmv, :mkl_dcscsv, :mkl_dcscmm, :mkl_dcscsm, :Float64),
(:mkl_ccscmv, :mkl_ccscsv, :mkl_ccscmm, :mkl_ccscsm, :Complex64),
(:mkl_zcscmv, :mkl_zcscsv, :mkl_zcscmm, :mkl_zcscsm, :Complex128))
(:mkl_ccscmv, :mkl_ccscsv, :mkl_ccscmm, :mkl_ccscsm, :ComplexF32),
(:mkl_zcscmv, :mkl_zcscsv, :mkl_zcscmm, :mkl_zcscsm, :ComplexF64))
@eval begin
function cscmv!(transa::Char, α::$T, matdescra::String,
A::SparseMatrixCSC{$T, BlasInt}, x::StridedVector{$T},
β::$T, y::StridedVector{$T})
_check_transa(transa)
_check_mat_mult_matvec(y, A, x, transa)
__counter[] += 1
ccall(($(string(mv)), :libmkl_rt), Void,
(Ptr{UInt8}, Ptr{BlasInt}, Ptr{BlasInt}, Ptr{$T},
ccall(($(string(mv)), :libmkl_rt), Cvoid,
(Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt}, Ref{$T},
Ptr{UInt8}, Ptr{$T}, Ptr{BlasInt}, Ptr{BlasInt},
Ptr{BlasInt}, Ptr{$T}, Ptr{$T}, Ptr{$T}),
&transa, &A.m, &A.n, &α,
matdescra, A.nzval, A.rowval, pointer(A.colptr, 1),
pointer(A.colptr, 2), x, &β, y)
Ptr{BlasInt}, Ptr{$T}, Ref{$T}, Ptr{$T}),
transa, A.m, A.n, α,
matdescra, A.nzval, A.rowval, A.colptr,
pointer(A.colptr, 2), x, β, y)
return y
end

Expand All @@ -56,15 +56,15 @@ function cscmm!(transa::Char, α::$T, matdescra::String,
mB, nB = size(B)
mC, nC = size(C)
__counter[] += 1
ccall(($(string(mm)), :libmkl_rt), Void,
(Ptr{UInt8}, Ptr{BlasInt}, Ptr{BlasInt}, Ptr{BlasInt},
Ptr{$T}, Ptr{UInt8}, Ptr{$T}, Ptr{BlasInt},
Ptr{BlasInt}, Ptr{BlasInt}, Ptr{$T}, Ptr{BlasInt},
Ptr{$T}, Ptr{$T}, Ptr{BlasInt}),
&transa, &A.m, &nC, &A.n,
&α, matdescra, A.nzval, A.rowval,
pointer(A.colptr, 1), pointer(A.colptr, 2), B, &mB,
&β, C, &mC)
ccall(($(string(mm)), :libmkl_rt), Cvoid,
(Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt}, Ref{BlasInt},
Ref{$T}, Ptr{UInt8}, Ptr{$T}, Ptr{BlasInt},
Ptr{BlasInt}, Ptr{BlasInt}, Ptr{$T}, Ref{BlasInt},
Ref{$T}, Ptr{$T}, Ref{BlasInt}),
transa, A.m, nC, A.n,
α, matdescra, A.nzval, A.rowval,
A.colptr, pointer(A.colptr, 2), B, mB,
β, C, mC)
return C
end

Expand All @@ -75,12 +75,12 @@ function cscsv!(transa::Char, α::$T, matdescra::String,
_check_transa(transa)
_check_mat_mult_matvec(y, A, x, transa)
__counter[] += 1
ccall(($(string(sv)), :libmkl_rt), Void,
(Ptr{UInt8}, Ptr{BlasInt}, Ptr{$T}, Ptr{UInt8},
ccall(($(string(sv)), :libmkl_rt), Cvoid,
(Ref{UInt8}, Ref{BlasInt}, Ref{$T}, Ptr{UInt8},
Ptr{$T}, Ptr{BlasInt}, Ptr{BlasInt}, Ptr{BlasInt},
Ptr{$T}, Ptr{$T}),
&transa, &A.m, &α, matdescra,
A.nzval, A.rowval, pointer(A.colptr, 1), pointer(A.colptr, 2),
transa, A.m, α, matdescra,
A.nzval, A.rowval, A.colptr, pointer(A.colptr, 2),
x, y)
return y
end
Expand All @@ -94,15 +94,15 @@ function cscsm!(transa::Char, α::$T, matdescra::String,
_check_transa(transa)
_check_mat_mult_matvec(C, A, B, transa)
__counter[] += 1
ccall(($(string(sm)), :libmkl_rt), Void,
(Ptr{UInt8}, Ptr{BlasInt}, Ptr{BlasInt}, Ptr{$T},
ccall(($(string(sm)), :libmkl_rt), Cvoid,
(Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt}, Ref{$T},
Ptr{UInt8}, Ptr{$T}, Ptr{BlasInt}, Ptr{BlasInt},
Ptr{BlasInt}, Ptr{$T}, Ptr{BlasInt}, Ptr{$T},
Ptr{BlasInt}),
&transa, &A.n, &nC, &α,
matdescra, A.nzval, A.rowval, pointer(A.colptr, 1),
pointer(A.colptr, 2), B, &mB, C,
&mC)
Ptr{BlasInt}, Ptr{$T}, Ref{BlasInt}, Ptr{$T},
Ref{BlasInt}),
transa, A.n, nC, α,
matdescra, A.nzval, A.rowval, A.colptr,
pointer(A.colptr, 2), B, mB, C,
mC)
return C
end

Expand Down
121 changes: 72 additions & 49 deletions src/BLAS/level_2_3/matmul.jl
Original file line number Diff line number Diff line change
@@ -1,77 +1,100 @@
import Base: *, A_mul_B!, At_mul_B!, Ac_mul_B!, Ac_mul_B, At_mul_B
import Base: \, A_ldiv_B!, At_ldiv_B!, Ac_ldiv_B!, Ac_ldiv_B, At_ldiv_B
import Base: *, \
import LinearAlgebra: mul!, ldiv!

_get_data(A::LowerTriangular) = tril(A.data)
_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

@compat 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}}}
_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 mat in (:StridedVector, :StridedMatrix)
for (trans, F!, F) in (('N', :A_mul_B! , :*),
('C', :Ac_mul_B!, :Ac_mul_B),
('T', :At_mul_B!, :At_mul_B))
for (tchar, ttype) in (('N', :()),
('C', :Adjoint),
('T', :Transpose))
AT = tchar == 'N' ? :(SparseMatrixCSC{$T,BlasInt}) : :($ttype{$T,SparseMatrixCSC{$T,BlasInt}})
@eval begin
function $F!::$T, A::SparseMatrixCSC{$T, BlasInt},
B::$mat{$T}, β::$T, C::$mat{$T})
isa(B,AbstractVector) ?
cscmv!($trans, α, matdescra(A), A, B, β, C) :
cscmm!($trans, α, matdescra(A), A, B, β, C)
function mul!::$T, adjA::$AT,
B::$mat{$T}, β::$T, C::$mat{$T})
A = _unwrap_adj(adjA)
if isa(B, AbstractVector)
return cscmv!($tchar, α, matdescra(A), A, B, β, C)
else
return cscmm!($tchar, α, matdescra(A), A, B, β, C)
end
end

function $F!(C::$mat{$T}, A::SparseMatrices{$T},
B::$mat{$T})
$F!(one($T), A, B, zero($T), C)
end
mul!(C::$mat{$T}, adjA::$AT, B::$mat{$T}) = mul!(one($T), adjA, B, zero($T), C)

function $F(A::SparseMatrices{$T}, B::$mat{$T})
isa(B,AbstractVector) ?
$F!(zeros($T, mkl_size($trans, A)[1]), A, B) :
$F!(zeros($T, mkl_size($trans, A)[1], size(B,2)), A, B)
function (*)(adjA::$AT, B::$mat{$T})
A = _unwrap_adj(adjA)
if isa(B,AbstractVector)
return mul!(zeros($T, mkl_size($tchar, A)[1]), adjA, B)
else
return mul!(zeros($T, mkl_size($tchar, A)[1], size(B,2)), adjA, B)
end
end
end

for w in (:Symmetric, :LowerTriangular, :UnitLowerTriangular, :UpperTriangular, :UnitUpperTriangular)
AT = tchar == 'N' ?
:($w{$T,SparseMatrixCSC{$T,BlasInt}}) :
:($ttype{$T,$w{$T,SparseMatrixCSC{$T,BlasInt}}})
@eval begin
function $F!::$T, A::$w{$T, SparseMatrixCSC{$T, BlasInt}},
function mul!::$T, adjA::$AT,
B::$mat{$T}, β::$T, C::$mat{$T})
isa(B,AbstractVector) ?
cscmv!($trans, α, matdescra(A), _get_data(A), B, β, C) :
cscmm!($trans, α, matdescra(A), _get_data(A), B, β, C)
A = _unwrap_adj(adjA)
if isa(B,AbstractVector)
return cscmv!($tchar, α, matdescra(A), _get_data(A), B, β, C)
else
return cscmm!($tchar, α, matdescra(A), _get_data(A), B, β, C)
end
end
end
end
end

for (trans, F!, F) in (('N', :A_ldiv_B! , :\),
('C', :Ac_ldiv_B!, :Ac_ldiv_B),
('T', :At_ldiv_B!, :At_ldiv_B))
for w in (:LowerTriangular, :UnitLowerTriangular, :UpperTriangular, :UnitUpperTriangular)
@eval begin
function $F!::$T, A::$w{$T, SparseMatrixCSC{$T, BlasInt}},
B::$mat{$T}, C::$mat{$T})
isa(B,AbstractVector) ?
cscsv!($trans, α, matdescra(A), _get_data(A), B, C) :
cscsm!($trans, α, matdescra(A), _get_data(A), B, C)
end
mul!(C::$mat{$T}, adjA::$AT, B::$mat{$T}) = mul!(one($T), adjA, B, zero($T), C)

function $F!(C::$mat{$T}, A::$w{$T, SparseMatrixCSC{$T, BlasInt}},
B::$mat{$T})
$F!(one($T), A, B, C)
function (*)(adjA::$AT, B::$mat{$T})
A = _unwrap_adj(adjA)
if isa(B,AbstractVector)
return mul!(zeros($T, mkl_size($tchar, A)[1]), adjA, B)
else
return mul!(zeros($T, mkl_size($tchar, A)[1], size(B,2)), adjA, B)
end
end
end

if w != :Symmetric
@eval begin
function ldiv!::$T, adjA::$AT,
B::$mat{$T}, C::$mat{$T})
A = _unwrap_adj(adjA)
if isa(B,AbstractVector)
return cscsv!($tchar, α, matdescra(A), _get_data(A), B, C)
else
return cscsm!($tchar, α, matdescra(A), _get_data(A), B, C)
end
end

ldiv!(C::$mat{$T}, A::$AT, B::$mat{$T}) =
ldiv!(one($T), A, B, C)

function $F(A::$w{$T, SparseMatrixCSC{$T, BlasInt}}, B::$mat{$T})
isa(B,AbstractVector) ?
$F!(zeros($T, size(A,1)), A, B) :
$F!(zeros($T, size(A,1), size(B,2)), A, B)
function (\)(A::$AT, B::$mat{$T})
if isa(B,AbstractVector)
return ldiv!(zeros($T, size(A,1)), A, B)
else
return ldiv!(zeros($T, size(A,1), size(B,2)), A, B)
end
end
end
end
end
Expand Down
2 changes: 0 additions & 2 deletions src/MKLSparse.jl
Original file line number Diff line number Diff line change
@@ -1,5 +1,3 @@
__precompile__()

module MKLSparse

function __init__()
Expand Down
47 changes: 22 additions & 25 deletions test/test_BLAS.jl
Original file line number Diff line number Diff line change
@@ -1,30 +1,27 @@
module Matdescra

using MKLSparse
using Base.Test
using Test, SparseArrays, LinearAlgebra

sA = sprand(5, 5, 0.01)
sS = sA'sA
sTl = tril(sS)
sTu = triu(sS)

@test MKLSparse.BLAS.matdescra(Base.LinAlg.Symmetric(sTl,:L)) == "SLNF"
@test MKLSparse.BLAS.matdescra(Base.LinAlg.Symmetric(sTu,:U)) == "SUNF"
@test MKLSparse.BLAS.matdescra(Base.LinAlg.LowerTriangular(sTl)) == "TLNF"
@test MKLSparse.BLAS.matdescra(Base.LinAlg.UpperTriangular(sTu)) == "TUNF"
@test MKLSparse.BLAS.matdescra(Base.LinAlg.UnitLowerTriangular(sTl)) == "TLUF"
@test MKLSparse.BLAS.matdescra(Base.LinAlg.UnitUpperTriangular(sTu)) == "TUUF"
@test MKLSparse.BLAS.matdescra(Symmetric(sTl,:L)) == "SLNF"
@test MKLSparse.BLAS.matdescra(Symmetric(sTu,:U)) == "SUNF"
@test MKLSparse.BLAS.matdescra(LowerTriangular(sTl)) == "TLNF"
@test MKLSparse.BLAS.matdescra(UpperTriangular(sTu)) == "TUNF"
@test MKLSparse.BLAS.matdescra(UnitLowerTriangular(sTl)) == "TLUF"
@test MKLSparse.BLAS.matdescra(UnitUpperTriangular(sTu)) == "TUUF"
@test MKLSparse.BLAS.matdescra(sA) == "GUUF"

end # module

module BLAS

import Base.LinAlg: UnitLowerTriangular, UnitUpperTriangular


using MKLSparse
using Base.Test
using Test, LinearAlgebra, SparseArrays

macro test_blas(ex)
return quote
Expand All @@ -43,32 +40,32 @@ end
@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.(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.(a.'*b - Array(a)'*b)) < 100*eps()
@test_blas maximum(abs.(transpose(a)*b - Array(a)'*b)) < 100*eps()
end
end

#?
@testset "complex matrix-vector multiplication" begin
for i = 1:5
a = speye(5) + im * 0.1*sprandn(5, 5, 0.2)
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(Complex128)
β = rand(Complex128)
α = 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.(a.'*b - Array(a).'*b)) < 100*eps())
@test_blas (maximum(abs.(A_mul_B!(similar(b), a, b) - Array(a)*b)) < 100*eps())
@test_blas (maximum(abs.(A_mul_B!(similar(c), a, c) - Array(a)*c)) < 100*eps())
@test_blas (maximum(abs.(At_mul_B!(similar(b), a, b) - Array(a).'*b)) < 100*eps())
@test_blas (maximum(abs.(At_mul_B!(similar(c), a, c) - Array(a).'*c)) < 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())

c = randn(6) + im*randn(6)
@test_throws DimensionMismatch a.'*c
@test_throws DimensionMismatch transpose(a)*c
@test_throws DimensionMismatch a.*c
@test_throws DimensionMismatch a.*c
end
Expand All @@ -78,11 +75,11 @@ end
n = 100
A = sprandn(n, n, 0.5) + sqrt(n)*I
b = rand(n)
symA = A + A.'
symA = A + transpose(A)
trilA = tril(A)
triuA = triu(A)
trilUA = tril(A, -1) + speye(A)
triuUA = triu(A, 1) + speye(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
Expand Down

0 comments on commit 3957a66

Please sign in to comment.