use libblastrampoline instead of libblas/liblapack (#157)

* don't import libblas/liblapack

* bump version to v0.6.9

* add tests

* tests for hesseneigvals

* tests for gemm

Author: jishnub

Project.toml

@@ -1,6 +1,6 @@
 name = "ApproxFunBase"
 uuid = "fbd15aa5-315a-5a7d-a8a4-24992e37be05"
-version = "0.6.8"
+version = "0.6.9"
 
 [deps]
 AbstractFFTs = "621f4979-c628-5d54-868e-fcf4e3e8185c"

src/LinearAlgebra/blas.jl

@@ -1,5 +1,7 @@
-import LinearAlgebra.BLAS: @blasfunc, libblas
-import LinearAlgebra.LAPACK: liblapack
+import LinearAlgebra.BLAS: @blasfunc
+const libblas = VERSION < v"1.8" ? LinearAlgebra.BLAS.libblas : "libblastrampoline"
+const liblapack = VERSION < v"1.8" ? LinearAlgebra.LAPACK.liblapack : "libblastrampoline"
+
 # Level 2
 ## mv
 ### gemv

src/LinearAlgebra/hesseneigs.jl

@@ -22,7 +22,7 @@ for (hseqr,elty) in ((:zhseqr_,:ComplexF64),)
 
         Ec='E'
         Nc='N'
-        ccall((@blasfunc($hseqr),LAPACK.liblapack),
+        ccall((@blasfunc($hseqr),liblapack),
             Nothing,
             (Ref{UInt8}, Ref{UInt8},
         Ref{BlasInt}, Ref{BlasInt}, Ref{BlasInt}, Ptr{$elty}, #A
@@ -59,7 +59,7 @@ for (hseqr,elty) in ((:dhseqr_,:Float64),)
         Ec='E'
         Nc='N'
         for i=1:2
-            ccall((@blasfunc($hseqr),LAPACK.liblapack),
+            ccall((@blasfunc($hseqr),liblapack),
                 Nothing,
                 (Ref{UInt8},Ref{UInt8},
             Ref{BlasInt},Ref{BlasInt},Ref{BlasInt},Ptr{$elty}, #A

test/runtests.jl

@@ -230,5 +230,54 @@ end
     @test a[2] ≈ exp(1)
 end
 
+@testset "BLAS/LAPACK" begin
+    @testset "gemv" begin
+        # test for the pointer versions, and assert that libblastrampoline works
+        @testset for T in [Float32, Float64, ComplexF32, ComplexF64]
+            a = zeros(T, 4)
+            b = zeros(T, 4)
+            A = Matrix{T}(I,4,4)
+            x = T[1:4;]
+            α, β = T(1.0), T(0.0)
+            ApproxFunBase.gemv!('N', α, A, x, β, a)
+            LinearAlgebra.BLAS.gemv!('N', α, A, x, β, b)
+            @test a == b == x
+            β = T(1.0)
+            ApproxFunBase.gemv!('N', α, A, x, β, a)
+            LinearAlgebra.BLAS.gemv!('N', α, A, x, β, b)
+            @test a == b == 2x
+        end
+    end
+
+    @testset "gemm" begin
+        # test for the pointer versions, and assert that libblastrampoline works
+        @testset for T in [Float32, Float64, ComplexF32, ComplexF64]
+            C1 = zeros(T, 4, 4)
+            C2 = zeros(T, 4, 4)
+            A = Matrix{T}(I,4,4)
+            B = reshape(T[1:16;], 4, 4)
+            α, β = T(1.0), T(0.0)
+            ApproxFunBase.gemm!('N', 'N', α, A, B, β, C1)
+            LinearAlgebra.BLAS.gemm!('N', 'N', α, A, B, β, C2)
+            @test C1 == C2 == B
+            β = T(1.0)
+            ApproxFunBase.gemm!('N', 'N', α, A, B, β, C1)
+            LinearAlgebra.BLAS.gemm!('N', 'N', α, A, B, β, C2)
+            @test C1 == C2 == 2B
+        end
+    end
+
+    @testset "hesseneigs" begin
+        A = Float64[1 4 2 3; 1 4 1 7; 0 2 3 4; 0 0 1 3]
+        λ1 = sort(ApproxFunBase.hesseneigvals(A), by = x->(real(x), imag(x)))
+        λ2 = eigvals(A)
+        @test λ1 ≈ λ2
+        B = ComplexF64.(A)
+        λ1 = sort(ApproxFunBase.hesseneigvals(B), by = x->(real(x), imag(x)))
+        λ2 = eigvals(B)
+        @test λ1 ≈ λ2
+    end
+end
+
 @time include("ETDRK4Test.jl")
 include("show.jl")