Update blas.jl

Update blas.jl 1. test negative stride more thoroughly 2. add missing test to `trmv` and `trsv` and put them together 3. add missing test to `symv` and `hemv`
JuliaLang · Nov 7, 2021 · 7ea2135 · 7ea2135
1 parent ad48bf1
commit 7ea2135
Showing 1 changed file with 130 additions and 125 deletions.
diff --git a/stdlib/LinearAlgebra/test/blas.jl b/stdlib/LinearAlgebra/test/blas.jl
@@ -88,50 +88,38 @@ end
             x2 = randn(elty, n)
             α  = rand(elty)
             β  = rand(elty)
-            @test BLAS.axpy!(α,copy(x1),copy(x2)) ≈ α*x1 + x2
-            @test BLAS.axpby!(α,copy(x1),β,copy(x2)) ≈ α*x1 + β*x2
+            for X1 in (x1, view(x1,n:-1:1)), X2 in (x2, view(x2, n:-1:1))
+                @test BLAS.axpy!(α,deepcopy(X1),deepcopy(X2)) ≈ α*X1 + X2
+                @test BLAS.axpby!(α,deepcopy(X1),β,deepcopy(X2)) ≈ α*X1 + β*X2
+            end
+            for ind1 in (1:n, n:-1:1), ind2 in (1:n, n:-1:1)
+                @test BLAS.axpy!(α,copy(x1),ind1,copy(x2),ind2) ≈ x2 + α*(ind1 == ind2 ? x1 : reverse(x1))
+            end
             @test_throws DimensionMismatch BLAS.axpy!(α, copy(x1), rand(elty, n + 1))
             @test_throws DimensionMismatch BLAS.axpby!(α, copy(x1), β, rand(elty, n + 1))
             @test_throws DimensionMismatch BLAS.axpy!(α, copy(x1), 1:div(n,2), copy(x2), 1:n)
             @test_throws ArgumentError BLAS.axpy!(α, copy(x1), 0:div(n,2), copy(x2), 1:(div(n, 2) + 1))
             @test_throws ArgumentError BLAS.axpy!(α, copy(x1), 1:div(n,2), copy(x2), 0:(div(n, 2) - 1))
-            @test BLAS.axpy!(α,copy(x1),1:n,copy(x2),1:n) ≈ x2 + α*x1
-            # negative stride test
-            @test BLAS.axpy!(α,copy(x1),n:-1:1,copy(x2),n:-1:1) ≈ x2 + α*x1
-            x1′, x2′ = @views x1[end:-1:1], x2[end:-1:1]
-            @test BLAS.axpy!(α, deepcopy(x1′), deepcopy(x2′)) ≈ α*x1′ + x2′
-            @test BLAS.axpy!(α, deepcopy(x1), deepcopy(x2′)) ≈ α*x1 + x2′
-            @test BLAS.axpy!(α, deepcopy(x1′), deepcopy(x2)) ≈ α*x1′ + x2
-            @test BLAS.axpby!(α, deepcopy(x1′), β, deepcopy(x2)) ≈ α*x1′ + β*x2
         end
         @testset "nrm2, iamax, and asum for StridedVectors" begin
             a = rand(elty,n)
-            b = view(a,2:2:n,1)
-            @test BLAS.nrm2(b) ≈ norm(b)
-            @test BLAS.asum(b) ≈ sum(fabs, b)
-            @test BLAS.iamax(b) == findmax(fabs, b)[2]
-            # negative stride test
-            c = view(a,n:-2:2)
-            @test BLAS.nrm2(c) ≈ norm(c)
-            @test BLAS.asum(c) ≈ sum(fabs, c)
-            @test BLAS.iamax(c) == 0
+            for ind in (2:2:n, n:-2:2)
+                b = view(a, ind, 1)
+                @test BLAS.nrm2(b) ≈ sqrt(sum(abs2, b))
+                @test BLAS.asum(b) ≈ sum(fabs, b)
+                @test BLAS.iamax(b) == findmax(fabs, b)[2] * (step(ind) >= 0)
+            end
         end
-        # scal
-        α = rand(elty)
-        a = rand(elty,n)
-        @test BLAS.scal(n,α,a,1) ≈ α * a
-        # negative stride test
-        @test BLAS.scal!(α, view(copy(a), n:-1:1)) ≈ α * reverse(a)
-
-        @testset "trsv" begin
-            A = triu(rand(elty,n,n))
-            @testset "Vector and SubVector" for x in (rand(elty, n), view(rand(elty,2n),1:2:2n), view(rand(elty,n),n:-1:1))
-                @test A\x ≈ BLAS.trsv('U','N','N', A, x) ≈ BLAS.trsv!('U','N','N', A, deepcopy(x))
-                @test_throws DimensionMismatch BLAS.trsv('U','N','N',A,Vector{elty}(undef,n+1))
+        @testset "scal" begin
+            α = rand(elty)
+            a = rand(elty,n)
+            @test BLAS.scal(n,α,a,1) ≈ α * a
+            for v in (a, view(a, n:-1:1))
+                @test BLAS.scal!(α, deepcopy(v)) ≈ α * v
             end
         end
-        @testset "ger, her, syr" for x in (rand(elty, n), view(rand(elty,2n), 1:2:2n), view(rand(elty,2n), 2n:-2:1)),
-            y in (rand(elty,n), view(rand(elty,3n), 1:3:3n), view(rand(elty,3n), 3n:-3:1))
+        @testset "ger, her, syr" for x in (rand(elty, n), view(rand(elty,2n), 1:2:2n), view(rand(elty,n), n:-1:1)),
+            y in (rand(elty,n), view(rand(elty,3n), 1:3:3n), view(rand(elty,2n), 2n:-2:2))
 
             A = rand(elty,n,n)
             α = rand(elty)
@@ -156,35 +144,63 @@ end
         @testset "copy" begin
             x1 = randn(elty, n)
             x2 = randn(elty, n)
-            @test x2 === BLAS.copyto!(x2, 1:n, x1, 1:n) == x1
+            for ind1 in (1:n, n:-1:1), ind2 in (1:n, n:-1:1)
+                @test x2 === BLAS.copyto!(x2, ind1, x1, ind2) == (ind1 == ind2 ? x1 : reverse(x1))
+            end
             @test_throws DimensionMismatch BLAS.copyto!(x2, 1:n, x1, 1:(n - 1))
             @test_throws ArgumentError BLAS.copyto!(x1, 0:div(n, 2), x2, 1:(div(n, 2) + 1))
             @test_throws ArgumentError BLAS.copyto!(x1, 1:(div(n, 2) + 1), x2, 0:div(n, 2))
-            @test x2 === BLAS.copyto!(x2, 1:n, x1, n:-1:1) == reverse(x1)
         end
-        # trmv
-        A = triu(rand(elty,n,n))
-        x = rand(elty,n)
-        @test BLAS.trmv('U','N','N',A,x) ≈ A*x
-        @test BLAS.trmv!('U','N','N',A,view(copy(x),n:-1:1)) ≈ A*x[n:-1:1]
+        @testset "trmv and trsv" begin
+            A = rand(elty,n,n)
+            x = rand(elty,n)
+            xerr = Vector{elty}(undef,n+1)
+            for uplo in ('U', 'L'), diag in ('U','N'), trans in ('N', 'T', 'C')
+                Wrapper = if uplo == 'U'
+                    diag == 'U' ? UnitUpperTriangular : UpperTriangular
+                else
+                    diag == 'U' ? UnitLowerTriangular : LowerTriangular
+                end
+                fun = trans == 'N' ? identity : trans == 'T' ? transpose : adjoint
+                fullA = collect(fun(Wrapper(A)))
+                @testset "trmv" begin
+                    @test BLAS.trmv(uplo,trans,diag,A,x) ≈ fullA * x
+                    @test_throws DimensionMismatch BLAS.trmv(uplo,trans,diag,A,xerr)
+                    for xx in (x, view(x, n:-1:1))
+                        @test BLAS.trmv!(uplo,trans,diag,A,deepcopy(xx)) ≈ fullA * xx
+                    end
+                end
+                @testset "trsv" begin
+                    @test BLAS.trsv(uplo,trans,diag,A,x) ≈ fullA \ x
+                    @test_throws DimensionMismatch BLAS.trsv(uplo,trans,diag,A,xerr)
+                    for xx in (x, view(x, n:-1:1))
+                        @test BLAS.trsv!(uplo,trans,diag,A,deepcopy(xx)) ≈ fullA \ xx
+                    end
+                end
+            end
+        end
         @testset "symmetric/Hermitian multiplication" begin
             x = rand(elty,n)
             A = rand(elty,n,n)
+            y = rand(elty, n)
+            α = randn(elty)
+            β = randn(elty)
             Aherm = A + A'
             Asymm = A + transpose(A)
-            @testset "symv and hemv" begin
-                @test BLAS.symv('U',Asymm,x) ≈ Asymm*x
-                @test BLAS.symv('U',Asymm,view(x,n:-1:1)) ≈ Asymm*x[n:-1:1]
-                @test BLAS.symv!('U',one(elty),Asymm,x,zero(elty),view(similar(x),n:-1:1)) ≈ Asymm*x
-                offsizevec, offsizemat = Array{elty}.(undef,(n+1, (n,n+1)))
-                @test_throws DimensionMismatch BLAS.symv!('U',one(elty),Asymm,x,one(elty),offsizevec)
-                @test_throws DimensionMismatch BLAS.symv('U',offsizemat,x)
+            offsizevec, offsizemat = Array{elty}.(undef,(n+1, (n,n+1)))
+            @testset "symv and hemv" for uplo in ('U', 'L')
+                @test BLAS.symv(uplo,Asymm,x) ≈ Asymm*x
+                for xx in (x, view(x, n:-1:1)), yy in (y, view(y, n:-1:1))
+                    @test BLAS.symv!(uplo,α,Asymm,xx,β,deepcopy(yy)) ≈ α * Asymm * xx + β * yy
+                end
+                @test_throws DimensionMismatch BLAS.symv!(uplo,α,Asymm,x,β,offsizevec)
+                @test_throws DimensionMismatch BLAS.symv(uplo,offsizemat,x)
                 if elty <: BlasComplex
-                    @test BLAS.hemv('U',Aherm,x) ≈ Aherm*x
-                    @test BLAS.hemv('U',Aherm,view(x,n:-1:1)) ≈ Aherm*x[n:-1:1]
-                    @test BLAS.hemv!('U',one(elty),Aherm,x,zero(elty),view(similar(x),n:-1:1)) ≈ Aherm*x
-                    @test_throws DimensionMismatch BLAS.hemv('U',offsizemat,x)
-                    @test_throws DimensionMismatch BLAS.hemv!('U',one(elty),Aherm,x,one(elty),offsizevec)
+                    @test BLAS.hemv(uplo,Aherm,x) ≈ Aherm*x
+                    for xx in (x, view(x, n:-1:1)), yy in (y, view(y, n:-1:1))
+                        @test BLAS.hemv!(uplo,α,Aherm,xx,β,deepcopy(yy)) ≈ α * Aherm * xx + β * yy                    end
+                    @test_throws DimensionMismatch BLAS.hemv(uplo,offsizemat,x)
+                    @test_throws DimensionMismatch BLAS.hemv!(uplo,one(elty),Aherm,x,one(elty),offsizevec)
                 end
             end
 
@@ -214,37 +230,27 @@ end
                 # Both matrix dimensions n coincide, as we have Hermitian matrices.
                 # Define the inputs and outputs of hpmv!, y = α*A*x+β*y
                 α = rand(elty)
-                M = rand(elty, n, n)
-                AL = Hermitian(M, :L)
-                AU = Hermitian(M, :U)
+                A = rand(elty, n, n)
                 x = rand(elty, n)
                 β = rand(elty)
                 y = rand(elty, n)
-
-                # Create lower triangular packing of AL
-                AP = elty[]
-                for j in 1:n, i in j:n
-                    push!(AP, AL[i,j])
+                for uplo in (:L, :U)
+                    Cuplo = String(uplo)[1]
+                    AH = Hermitian(A, uplo)
+                    # Create lower/upper triangular packing of AL
+                    AP = elty[]
+                    for j in 1:n, i in (uplo == :L ? (j:n) : (1:j))
+                        push!(AP, AH[i,j])
+                    end
+                    for xx in (x, view(x,n:-1:1)), yy in (y, view(y,n:-1:1))
+                        @test BLAS.hpmv!(Cuplo, α, AP, xx, β, deepcopy(yy)) ≈ α*AH*xx + β*yy
+                    end
+                    AP′ = view(zeros(elty, n*(n+1)),1:2:n*(n+1))
+                    @test_throws ArgumentError BLAS.hpmv!(Cuplo, α, AP′, x, β, y)
+                    AP′ = view(AP, 1:length(AP′) - 1)
+                    @test_throws DimensionMismatch BLAS.hpmv!(Cuplo, α, AP′, x, β, y)
+                    @test_throws DimensionMismatch BLAS.hpmv!(Cuplo, α, AP′, x, β, view(y,1:n-1))
                 end
-                @test BLAS.hpmv!('L', α, AP, x, β, copy(y)) ≈ α*AL*x + β*y
-                @test BLAS.hpmv!('L', α, AP, view(x,n:-1:1), β, copy(y)) ≈
-                    BLAS.hpmv!('L', α, AP, x[n:-1:1], β, copy(y))
-
-                # Create upper triangular packing of AU
-                AP = elty[]
-                for j in 1:n, i in 1:j
-                    push!(AP, AU[i,j])
-                end
-                @test BLAS.hpmv!('U', α, AP, x, β, copy(y)) ≈ α*AU*x + β*y
-                @test BLAS.hpmv!('U', α, AP, view(x,n:-1:1), β, copy(y)) ≈
-                    BLAS.hpmv!('U', α, AP, x[n:-1:1], β, copy(y))
-
-                # inputs check
-                AP′ = view(zeros(elty, n*(n+1)),1:2:n*(n+1))
-                @test_throws ArgumentError BLAS.hpmv!('L', α, AP′, x, β, y)
-                AP′ = view(AP, 1:length(AP′) - 1)
-                @test_throws DimensionMismatch BLAS.hpmv!('L', α, AP′, x, β, y)
-                @test_throws DimensionMismatch BLAS.hpmv!('L', α, AP′, x, β, view(y,1:n-1))
             end
         end
 
@@ -254,83 +260,82 @@ end
                 # Both matrix dimensions n coincide, as we have symmetric matrices.
                 # Define the inputs and outputs of spmv!, y = α*A*x+β*y
                 α = rand(elty)
-                M = rand(elty, n, n)
-                AL = Symmetric(M, :L)
-                AU = Symmetric(M, :U)
+                A = rand(elty, n, n)
                 x = rand(elty, n)
                 β = rand(elty)
                 y = rand(elty, n)
-
-                # Create lower triangular packing of AL
-                AP = elty[]
-                for j in 1:n, i in j:n
-                    push!(AP, AL[i,j])
+                for uplo in (:L, :U)
+                    Cuplo = String(uplo)[1]
+                    AS = Symmetric(A, uplo)
+                    # Create lower/upper triangular packing of AL
+                    AP = elty[]
+                    for j in 1:n, i in (uplo == :L ? (j:n) : (1:j))
+                        push!(AP, AS[i,j])
+                    end
+                    for xx in (x, view(x,n:-1:1)), yy in (y, view(y,n:-1:1))
+                        @test BLAS.spmv!(Cuplo, α, AP, xx, β, deepcopy(yy)) ≈ α*AS*xx + β*yy
+                    end
+                    AP′ = view(zeros(elty, n*(n+1)),1:2:n*(n+1))
+                    @test_throws ArgumentError BLAS.spmv!(Cuplo, α, AP′, x, β, y)
+                    AP′ = view(AP, 1:length(AP′) - 1)
+                    @test_throws DimensionMismatch BLAS.spmv!(Cuplo, α, AP′, x, β, y)
+                    @test_throws DimensionMismatch BLAS.spmv!(Cuplo, α, AP′, x, β, view(y,1:n-1))
                 end
-                @test BLAS.spmv!('L', α, AP, x, β, copy(y)) ≈ α*AL*x + β*y
-                @test BLAS.spmv!('L', α, AP, view(x,n:-1:1), β, copy(y)) ≈
-                    BLAS.spmv!('L', α, AP, x[n:-1:1], β, copy(y))
-
-                y_result_julia_upper = α*AU*x + β*y
-
-                # Create upper triangular packing of AU
-                AP = elty[]
-                for j in 1:n, i in 1:j
-                    push!(AP, AU[i,j])
-                end
-                @test BLAS.spmv!('U', α, AP, x, β, copy(y)) ≈ α*AU*x + β*y
-                @test BLAS.spmv!('U', α, AP, view(x,n:-1:1), β, copy(y)) ≈
-                    BLAS.spmv!('U', α, AP, x[n:-1:1], β, copy(y))
-
-                # inputs check
-                AP′ = view(zeros(elty, n*(n+1)),1:2:n*(n+1))
-                @test_throws ArgumentError BLAS.spmv!('L', α, AP′, x, β, y)
-                AP′ = view(AP, 1:length(AP′) - 1)
-                @test_throws DimensionMismatch BLAS.spmv!('L', α, AP′, x, β, y)
-                @test_throws DimensionMismatch BLAS.spmv!('L', α, AP′, x, β, view(y,1:n-1))
             end
         end
-
-        #trsm
-        A = triu(rand(elty,n,n))
-        B = rand(elty,(n,n))
-        @test BLAS.trsm('L','U','N','N',one(elty),A,B) ≈ A\B
-
+        @testset "trsm" begin
+            A = triu(rand(elty,n,n))
+            B = rand(elty,(n,n))
+            @test BLAS.trsm('L','U','N','N',one(elty),A,B) ≈ A\B
+        end
         #will work for SymTridiagonal,Tridiagonal,Bidiagonal!
         @testset "banded matrix mv" begin
             @testset "gbmv" begin
-                TD  = Tridiagonal(rand(elty,n-1),rand(elty,n),rand(elty,n-1))
-                x   = rand(elty,n)
+                TD = Tridiagonal(rand(elty,n-1),rand(elty,n),rand(elty,n-1))
+                x  = rand(elty, n)
                 #put TD into the BLAS format!
                 fTD = zeros(elty,3,n)
                 fTD[1,2:n] = TD.du
                 fTD[2,:] = TD.d
                 fTD[3,1:n-1] = TD.dl
                 @test BLAS.gbmv('N',n,1,1,fTD,x) ≈ TD*x
-                @test BLAS.gbmv('N',n,1,1,fTD,view(x,n:-1:1)) ≈ TD*x[n:-1:1]
-                @test BLAS.gbmv!('N',n,1,1,one(elty),fTD,x,zero(elty),view(similar(x),n:-1:1)) ≈ TD*x
+                y = rand(elty, n)
+                α = randn(elty)
+                β = randn(elty)
+                for xx in (x, view(x, n:-1:1)), yy in (y, view(y, n:-1:1))
+                    @test BLAS.gbmv!('N',n,1,1,α,fTD,xx,β,deepcopy(yy)) ≈ α * TD * xx + β * yy
+                end
             end
             #will work for SymTridiagonal only!
-            @testset "sbmv" begin
+            @testset "sbmv and hbmv" begin
+                x   = rand(elty,n)
                 if elty <: BlasReal
                     ST  = SymTridiagonal(rand(elty,n),rand(elty,n-1))
-                    x   = rand(elty,n)
                     #put TD into the BLAS format!
                     fST = zeros(elty,2,n)
                     fST[1,2:n] = ST.ev
                     fST[2,:] = ST.dv
                     @test BLAS.sbmv('U',1,fST,x) ≈ ST*x
-                    @test BLAS.sbmv('U',1,fST,view(x, n:-1:1)) ≈ ST*x[n:-1:1]
-                    @test BLAS.sbmv!('U',1, one(elty), fST, x, zero(elty), view(similar(x),n:-1:1)) ≈ ST*x
+                    y = rand(elty, n)
+                    α = randn(elty)
+                    β = randn(elty)
+                    for xx in (x, view(x, n:-1:1)), yy in (y, view(y, n:-1:1))
+                        @test BLAS.sbmv!('U',1,α,fST,xx,β,deepcopy(yy)) ≈ α * ST * xx + β * yy
+                    end
                 else
-                    dv = real(rand(elty,n))
+                    dv = rand(real(elty),n)
                     ev = rand(elty,n-1)
                     bH = zeros(elty,2,n)
                     bH[1,2:n] = ev
                     bH[2,:] = dv
                     fullH = diagm(0 => dv, -1 => conj(ev), 1 => ev)
                     @test BLAS.hbmv('U',1,bH,x) ≈ fullH*x
-                    @test BLAS.hbmv('U',1,bH,view(x,n:-1:1)) ≈ fullH*x[n:-1:1]
-                    @test BLAS.hbmv!('U',1, one(elty), bH, x, zero(elty), view(similar(x),n:-1:1)) ≈ fullH*x
+                    y = rand(elty, n)
+                    α = randn(elty)
+                    β = randn(elty)
+                    for xx in (x, view(x, n:-1:1)), yy in (y, view(y, n:-1:1))
+                        @test BLAS.hbmv!('U',1,α,bH,xx,β,deepcopy(yy)) ≈ α * fullH * xx + β * yy
+                    end
                 end
             end
         end