Updated for Julia 1.0, added (limited) auto detection of libraries and avx support

Project.toml

@@ -0,0 +1,10 @@
+name = "VML"
+uuid = "c8ce9da6-5d36-5c03-b118-5a70151be7bc"
+
+[deps]
+CpuId = "adafc99b-e345-5852-983c-f28acb93d879"
+Libdl = "8f399da3-3557-5675-b5ff-fb832c97cbdb"
+SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
+
+[compat]
+julia = "≥ 0.7 1.0"

README.md

@@ -1,25 +1,27 @@
 # VML
 
 This package provides bindings to the Intel Vector Math Library for
-arithmetic and transcendental functions. It is often substantially
-faster than using Julia's built-in functions.
+arithmetic and transcendental functions. Especially for large vectors it is often substantially faster than broadcasting Julia's built-in functions.
 
-## Using VML.jl
+## Setting up VML.jl
 
 To use VML.jl, you must have the Intel Vector Math Library installed.
-This is included in [MKL](http://software.intel.com/en-us/intel-mkl),
-which is free for non-commercial use. You must also copy/symlink the
-appropriate shared library to a directory known to the linker (e.g.
-`/usr/local/lib`) or you must modify the path to `lib` in `src/VML.jl`.
+For this you have two options. First, you can install (and build) [MKL.jl](https://github.com/JuliaComputing/MKL.jl), which will add the necessary libraries to your Julia install. This will change your Julia system image however, so if you would prefer not to do that you can get the stand-alone [MKL](http://software.intel.com/en-us/intel-mkl),
+which is free for non-commercial use. The default install location will be detected automatically (currently only `opt/intel/mkl/lib` on Unix). If you have chosen a different location or have multiple versions of MKL, you can also specify the environmental variable 
+```
+    ENV["MKL_SL"] = <.../lib>
+```
+The specified folder should contain `libmkl_rt`, `libmkl_core` and `libmkl_vml_avx` or `libmkl_vml_avx2`, with file endings appropriate for your operating system. 
+The definition of `ENV["MKL_SL"]` will take precedence even if `MKL.jl` or default standalone MKL are installed. 
+
+Using [CpuId.jl](https://github.com/m-j-w/CpuId.jl), VML.jl detects if your processor supports the newer `avx2` instructions, and if not default to `libmkl_vml_avx`. If your system does not have AVX this package will currently not work for you.
 
-Currently, VML.jl is configured to use `libmkl_vml_avx`, which requires
-AVX support. If your system does not have AVX (e.g., most pre-Sandy
-Bridge systems), you will need to modify the `const lib` declaration at
-the top of `src/VML.jl`. Future versions of VML.jl may automatically
-detect CPU architecture.
+If the CPU feature detection does not work for you, please open an issue. 
+
+## Using VML.jl
 
 After loading VML.jl, vector calls to functions listed below will
-automatically use VML instead of openlibm when possible.
+automatically use VML instead of openlibm when possible. Note that most function currently do not have a vectorized version (e.g. you call `sin.(rand(300))`), so there should be no conflict. If there is let me know. Updated and conflict tested exported functions are planned for the future.
 
 By default, VML uses `VML_HA` mode, which corresponds to an accuracy of
 <1 ulp, matching the accuracy of Julia's built-in openlibm
@@ -30,7 +32,7 @@ regarding these options is available on
 [Intel's website](http://software.intel.com/sites/products/documentation/hpc/mkl/vml/vmldata.htm).
 
 ## Performance
-
+(These results are currently outdated will be updated in due course)
 ![VML Performance Comparison](/benchmark/performance.png)
 
 ![VML Complex Performance Comparison](/benchmark/performance_complex.png)
@@ -50,7 +52,7 @@ these are not yet implemented in VML.jl.
 ### Unary functions
 
 Allocating forms have signature `f(A)`. Mutating forms have signatures
-`f!(A)` (in place) and `f!(out, A)` (out of place).
+`f!(A)` (in place) and `f!(out, A)` (out of place). The last 9 functions have been moved from Base to `SpecialFunctions.jl` or have no Base equivalent. 
 
 Allocating | Mutating
 -----------|---------
@@ -82,6 +84,8 @@ Allocating | Mutating
 `erfc`     | `erfc!`
 `erfinv`   | `erfinv!`
 `efcinv`   | `efcinv!`
+`gamma`    | `gamma!`
+`lgamma`   | `lgamma!`
 `inv_cbrt` | `inv_cbrt!`
 `inv_sqrt` | `inv_sqrt!`
 `pow2o3`   | `pow2o3!`
@@ -97,5 +101,9 @@ Allocating | Mutating
 -----------|---------
 `atan2`    | `atan2!`
 `hypot`    | `hypot!`
-`.^`       | `pow!`
-`./`       | `divide!`
+`pow`       | `pow!`
+`divide`       | `divide!`
+
+
+### Next steps
+Next steps for this package are proper Windows support, writing up proper testing to make sure each advertised function actually works as expected and build testing on travis.

deps/.gitignore

@@ -0,0 +1,2 @@
+deps.jl
+build.log

deps/build.jl

@@ -0,0 +1,81 @@
+
+using CpuId
+using Libdl
+
+
+function defMKLpath()
+    if Sys.isunix()
+        return "/opt/intel/mkl/lib"
+    # elseif Sys.iswindows()
+    #     return
+    else
+        return ""
+    end
+end
+
+
+
+### First Option, check for user preference
+if haskey(ENV, "MKL_SL") 
+    basepath = ENV["MKL_SL"]
+    if !any([occursin("libmkl_rt",x) for x in readdir(basepath)])
+        error(""" ENV["MKL_SL"] is set, but does not point to folder containing MKL files """)
+    end
+    println("""Found MKL via ENV["MKL_SL"]
+    basepath set to $basepath""")
+
+### Next, check if MKL.jl is already installed. 
+# Don't want to add to deps, so bit hacky
+else
+    foundMKLjl = [isdir(joinpath(testdir, "packages/MKL")) for testdir in DEPOT_PATH]
+
+    if any(foundMKLjl)
+        pkgPath = joinpath(DEPOT_PATH[findfirst(foundMKLjl)], "packages/MKL/")
+
+        shortID = readdir(pkgPath)[1]
+
+        basepath = joinpath(pkgPath, shortID, "deps/usr/lib")
+
+        if isdir(basepath)
+            println("Found MKL.jl package
+            basepath set to $basepath")
+        else
+            error("Found MKL package folder, but MKL libraries at deps/usr/lib are missing
+            Please rebuild MKL")
+        end
+
+    elseif isdir(defMKLpath())
+        basepath = defMKLpath()
+
+    else
+        error("Could not find VML shared libraries
+            Check github.com/.... for details on obtaining them")
+
+    end
+end
+
+
+if cpufeature(:AVX2)
+    lib = joinpath(basepath, "libmkl_vml_avx2")
+    println("AVX2 support detected, vml_avx2 selected")
+else
+    lib = joinpath(basepath, "libmkl_vml_avx")
+    println("AVX2 support missing, vml_avx selected")
+end
+
+rtlib = joinpath(basepath, "libmkl_rt")
+corelib = joinpath(basepath, "libmkl_core")
+
+
+depsjl_path = joinpath(@__DIR__, "deps.jl")
+open(depsjl_path, "w") do depsjl_file
+    println(depsjl_file, strip("""
+        ## This file was autogenerated by build.jl.
+        ## Do not edit.
+        import Libdl
+
+        const lib = "$lib"
+        const rtlib = "$rtlib"
+        const corelib = "$corelib"
+        """))
+end

src/VML.jl

@@ -1,110 +1,17 @@
-module VML
+__precompile__()
 
-import Base: .^, ./
+module VML
 
+# import Base: .^, ./
+using SpecialFunctions
+using Libdl
 # TODO detect CPU architecture
-const lib = :libmkl_vml_avx
-Libdl.dlopen(:libmkl_rt)
-
-immutable VMLAccuracy
-    mode::UInt
-end
-const VML_LA = VMLAccuracy(0x00000001)
-const VML_HA = VMLAccuracy(0x00000002)
-const VML_EP = VMLAccuracy(0x00000003)
-Base.show(io::IO, m::VMLAccuracy) = print(io, m == VML_LA ? "VML_LA" :
-                                              m == VML_HA ? "VML_HA" : "VML_EP")
-vml_get_mode() = ccall((:_vmlGetMode, lib), Cuint, ())
-vml_set_mode(mode::Integer) = (ccall((:_vmlSetMode, lib), Cuint, (UInt,), mode); nothing)
+# include("libdetect.jl")
+include(joinpath(dirname(@__DIR__), "deps/deps.jl"))
 
-vml_set_accuracy(m::VMLAccuracy) = vml_set_mode((vml_get_mode() & ~0x03) | m.mode)
-vml_get_accuracy() = VMLAccuracy(vml_get_mode() & 0x3)
-
-vml_set_mode((vml_get_mode() & ~0x0000FF00))
-function vml_check_error()
-    vml_error = ccall((:_vmlClearErrStatus, lib), Cint, ())
-    if vml_error != 0
-        if vml_error == 1
-            throw(DomainError())
-        elseif vml_error == 2 || vml_error == 3 || vml_error == 4
-            # Singularity, overflow, or underflow
-            # I don't think Base throws on these
-        elseif vml_error == 1000
-            warn("VML does not support $(vml_get_accuracy); lower accuracy used instead")
-        else
-            error("an unexpected error occurred in VML ($vml_error)")
-        end
-    end
-end
-
-function vml_prefix(t::DataType)
-    if t == Float32
-        return "_vmls"
-    elseif t == Float64
-        return "_vmld"
-    elseif t == Complex{Float32}
-        return "_vmlc"
-    elseif t == Complex{Float64}
-        return "_vmlz"
-    end
-    error("unknown type $t")
-end
+include("setup.jl")
 
-function def_unary_op(tin, tout, jlname, jlname!, mklname)
-    mklfn = Base.Meta.quot(Symbol("$(vml_prefix(tin))$mklname"))
-    exports = Symbol[]
-    isa(jlname, Expr) || push!(exports, jlname)
-    isa(jlname!, Expr) || push!(exports, jlname!)
-    @eval begin
-        $(isempty(exports) ? nothing : Expr(:export, exports...))
-        function $(jlname!){N}(out::Array{$tout,N}, A::Array{$tin,N})
-            size(out) == size(A) || throw(DimensionMismatch())
-            ccall(($mklfn, lib), Void, (Int, Ptr{$tin}, Ptr{$tout}), length(A), A, out)
-            vml_check_error()
-            out
-        end
-        $(if tin == tout
-            quote
-                function $(jlname!)(A::Array{$tin})
-                    ccall(($mklfn, lib), Void, (Int, Ptr{$tin}, Ptr{$tout}), length(A), A, A)
-                    vml_check_error()
-                    A
-                end
-            end
-        end)
-        function $(jlname)(A::Array{$tin})
-            out = similar(A, $tout)
-            ccall(($mklfn, lib), Void, (Int, Ptr{$tin}, Ptr{$tout}), length(A), A, out)
-            vml_check_error()
-            out
-        end
-    end
-end
-
-function def_binary_op(tin, tout, jlname, jlname!, mklname, broadcast)
-    mklfn = Base.Meta.quot(Symbol("$(vml_prefix(tin))$mklname"))
-    exports = Symbol[]
-    isa(jlname, Expr) || push!(exports, jlname)
-    isa(jlname!, Expr) || push!(exports, jlname!)
-    @eval begin
-        $(isempty(exports) ? nothing : Expr(:export, exports...))
-        function $(jlname!){N}(out::Array{$tout,N}, A::Array{$tin,N}, B::Array{$tin,N})
-            size(out) == size(A) == size(B) || $(broadcast ? :(return broadcast!($jlname, out, A, B)) : :(throw(DimensionMismatch())))
-            ccall(($mklfn, lib), Void, (Int, Ptr{$tin}, Ptr{$tin}, Ptr{$tout}), length(A), A, B, out)
-            vml_check_error()
-            out
-        end
-        function $(jlname){N}(A::Array{$tout,N}, B::Array{$tin,N})
-            size(A) == size(B) || $(broadcast ? :(return broadcast($jlname, A, B)) : :(throw(DimensionMismatch())))
-            out = similar(A)
-            ccall(($mklfn, lib), Void, (Int, Ptr{$tin}, Ptr{$tin}, Ptr{$tout}), length(A), A, B, out)
-            vml_check_error()
-            out
-        end
-    end
-end
-
-for t in (Float32, Float64, Complex64, Complex128)
+for t in (Float32, Float64, ComplexF32, ComplexF64)
     # Unary, real or complex
     def_unary_op(t, t, :(Base.acos), :acos!, :Acos)
     def_unary_op(t, t, :(Base.asin), :asin!, :Asin)
@@ -114,9 +21,9 @@ for t in (Float32, Float64, Complex64, Complex128)
     def_unary_op(t, t, :(Base.exp), :exp!, :Exp)
     def_unary_op(t, t, :(Base.log), :log!, :Ln)
 
-    # Binary, real or complex
-    def_binary_op(t, t, :(.^), :pow!, :Pow, true)
-    def_binary_op(t, t, :(./), :divide!, :Div, true)
+    # # Binary, real or complex
+    def_binary_op(t, t, :pow, :pow!, :Pow, true)
+    def_binary_op(t, t, :divide, :divide!, :Div, true)
 end
 
 for t in (Float32, Float64)
@@ -130,12 +37,14 @@ for t in (Float32, Float64)
     def_unary_op(t, t, :(Base.floor), :floor!, :Floor)
     def_unary_op(t, t, :(Base.round), :round!, :Round)
     def_unary_op(t, t, :(Base.trunc), :trunc!, :Trunc)
-    def_unary_op(t, t, :(Base.erf), :erf!, :Erf)
-    def_unary_op(t, t, :(Base.erfc), :erfc!, :Erfc)
-    def_unary_op(t, t, :(Base.erfinv), :erfinv!, :ErfInv)
-    def_unary_op(t, t, :(Base.erfcinv), :erfcinv!, :ErfcInv)
-    def_unary_op(t, t, :(Base.lgamma), :lgamma!, :LGamma)
-    def_unary_op(t, t, :(Base.gamma), :gamma!, :TGamma)
+
+    # now in SpecialFunctions (make smart, maybe?)
+    def_unary_op(t, t, :(SpecialFunctions.erf), :erf!, :Erf)
+    def_unary_op(t, t, :(SpecialFunctions.erfc), :erfc!, :Erfc)
+    def_unary_op(t, t, :(SpecialFunctions.erfinv), :erfinv!, :ErfInv)
+    def_unary_op(t, t, :(SpecialFunctions.erfcinv), :erfcinv!, :ErfcInv)
+    def_unary_op(t, t, :(SpecialFunctions.lgamma), :lgamma!, :LGamma)
+    def_unary_op(t, t, :(SpecialFunctions.gamma), :gamma!, :TGamma)
     # Not in Base
     def_unary_op(t, t, :inv_cbrt, :inv_cbrt!, :InvCbrt)
     def_unary_op(t, t, :inv_sqrt, :inv_sqrt!, :InvSqrt)
@@ -154,31 +63,31 @@ for t in (Float32, Float64)
     def_unary_op(t, t, :(Base.tanh), :tanh!, :Tanh)
     def_unary_op(t, t, :(Base.log10), :log10!, :Log10)
 
-    # .^ to scalar power
-    mklfn = Base.Meta.quot(Symbol("$(vml_prefix(t))Powx"))
-    @eval begin
-        export pow!
-        function pow!{N}(out::Array{$t,N}, A::Array{$t,N}, b::$t)
-            size(out) == size(A) || throw(DimensionMismatch())
-            ccall(($mklfn, lib), Void, (Int, Ptr{$t}, $t, Ptr{$t}), length(A), A, b, out)
-            vml_check_error()
-            out
-        end
-        function (.^){N}(A::Array{$t,N}, b::$t)
-            out = similar(A)
-            ccall(($mklfn, lib), Void, (Int, Ptr{$t}, $t, Ptr{$t}), length(A), A, b, out)
-            vml_check_error()
-            out
-        end
-    end
-
-    # Binary, real-only
-    def_binary_op(t, t, :(Base.atan2), :atan2!, :Atan2, false)
+    # # .^ to scalar power
+    # mklfn = Base.Meta.quot(Symbol("$(vml_prefix(t))Powx"))
+    # @eval begin
+    #     export pow!
+    #     function pow!{N}(out::Array{$t,N}, A::Array{$t,N}, b::$t)
+    #         size(out) == size(A) || throw(DimensionMismatch())
+    #         ccall(($mklfn, lib), Nothing, (Int, Ptr{$t}, $t, Ptr{$t}), length(A), A, b, out)
+    #         vml_check_error()
+    #         out
+    #     end
+    #     function (.^){N}(A::Array{$t,N}, b::$t)
+    #         out = similar(A)
+    #         ccall(($mklfn, lib), Nothing, (Int, Ptr{$t}, $t, Ptr{$t}), length(A), A, b, out)
+    #         vml_check_error()
+    #         out
+    #     end
+    # end
+
+    # # Binary, real-only
+    def_binary_op(t, t, :(Base.atan), :atan!, :Atan, false)
     def_binary_op(t, t, :(Base.hypot), :hypot!, :Hypot, false)
 
     # Unary, complex-only
     def_unary_op(t, Complex{t}, :(Base.cis), :cis!, :CIS)
-    # def_unary_op(Complex{t}, Complex{t}, :(Base.conj), :conj!, :Conj)
+    def_unary_op(Complex{t}, Complex{t}, :(Base.conj), :conj!, :Conj)
     def_unary_op(Complex{t}, t, :(Base.abs), :abs!, :Abs)
     def_unary_op(Complex{t}, t, :(Base.angle), :angle!, :Arg)
 
@@ -193,4 +102,5 @@ for t in (Float32, Float64)
 end
 
 export VML_LA, VML_HA, VML_EP, vml_set_accuracy, vml_get_accuracy
+
 end