diff --git a/stdlib/LinearAlgebra/src/blas.jl b/stdlib/LinearAlgebra/src/blas.jl index 661e9e2b15617..a0db0518e25ba 100644 --- a/stdlib/LinearAlgebra/src/blas.jl +++ b/stdlib/LinearAlgebra/src/blas.jl @@ -159,6 +159,18 @@ end # Level 1 +isdense(x) = x isa DenseArray +isdense(x::Base.FastContiguousSubArray) = isdense(parent(x)) +isdense(x::Base.ReshapedArray) = isdense(parent(x)) +isdense(x::Base.ReinterpretArray) = isdense(parent(x)) +@inline function ptrst1(x::AbstractArray) + isdense(x) && return pointer(x), 1 # simpify runtime check when possibe + ndims(x) == 1 || strides(x) == Base.size_to_strides(strides(x, 1), size(x)...) || + throw(ArgumentError("only support vector like inputs")) + st = stride(x, 1) + ptr = st >= 0 ? pointer(x) : pointer(x, lastindex(x)) + ptr, st +end ## copy """ @@ -249,7 +261,10 @@ for (fname, elty) in ((:dscal_,:Float64), DX end - scal!(DA::$elty, DX::AbstractArray{$elty}) = scal!(length(DX),DA,DX,stride(DX,1)) + scal!(DA::$elty, DX::AbstractArray{$elty}) = let (p, st) = ptrst1(DX) + GC.@preserve DX scal!(length(DX), DA, p, abs(st)) + DX + end end end scal(n, DA, DX, incx) = scal!(n, DA, copy(DX), incx) @@ -353,75 +368,18 @@ for (fname, elty) in ((:cblas_zdotu_sub,:ComplexF64), end end -@inline function _dot_length_check(x,y) - n = length(x) - if n != length(y) - throw(DimensionMismatch("dot product arguments have lengths $(length(x)) and $(length(y))")) - end - n -end - for (elty, f) in ((Float32, :dot), (Float64, :dot), (ComplexF32, :dotc), (ComplexF64, :dotc), (ComplexF32, :dotu), (ComplexF64, :dotu)) @eval begin - function $f(x::DenseArray{$elty}, y::DenseArray{$elty}) - n = _dot_length_check(x,y) - $f(n, x, 1, y, 1) - end - - function $f(x::StridedVector{$elty}, y::DenseArray{$elty}) - n = _dot_length_check(x,y) - xstride = stride(x,1) - ystride = stride(y,1) - x_delta = xstride < 0 ? n : 1 - GC.@preserve x $f(n,pointer(x,x_delta),xstride,y,ystride) - end - - function $f(x::DenseArray{$elty}, y::StridedVector{$elty}) - n = _dot_length_check(x,y) - xstride = stride(x,1) - ystride = stride(y,1) - y_delta = ystride < 0 ? n : 1 - GC.@preserve y $f(n,x,xstride,pointer(y,y_delta),ystride) - end - - function $f(x::StridedVector{$elty}, y::StridedVector{$elty}) - n = _dot_length_check(x,y) - xstride = stride(x,1) - ystride = stride(y,1) - x_delta = xstride < 0 ? n : 1 - y_delta = ystride < 0 ? n : 1 - GC.@preserve x y $f(n,pointer(x,x_delta),xstride,pointer(y,y_delta),ystride) + function $f(x::AbstractVector{$elty}, y::AbstractVector{$elty}) + n, m = length(x), length(y) + n == m || throw(DimensionMismatch("dot product arguments have lengths $n and $m")) + GC.@preserve x y $f(n, ptrst1(x)..., ptrst1(y)...) end end end -function dot(DX::Union{DenseArray{T},AbstractVector{T}}, DY::Union{DenseArray{T},AbstractVector{T}}) where T<:BlasReal - require_one_based_indexing(DX, DY) - n = length(DX) - if n != length(DY) - throw(DimensionMismatch("dot product arguments have lengths $(length(DX)) and $(length(DY))")) - end - return dot(n, DX, stride(DX, 1), DY, stride(DY, 1)) -end -function dotc(DX::Union{DenseArray{T},AbstractVector{T}}, DY::Union{DenseArray{T},AbstractVector{T}}) where T<:BlasComplex - require_one_based_indexing(DX, DY) - n = length(DX) - if n != length(DY) - throw(DimensionMismatch("dot product arguments have lengths $(length(DX)) and $(length(DY))")) - end - return dotc(n, DX, stride(DX, 1), DY, stride(DY, 1)) -end -function dotu(DX::Union{DenseArray{T},AbstractVector{T}}, DY::Union{DenseArray{T},AbstractVector{T}}) where T<:BlasComplex - require_one_based_indexing(DX, DY) - n = length(DX) - if n != length(DY) - throw(DimensionMismatch("dot product arguments have lengths $(length(DX)) and $(length(DY))")) - end - return dotu(n, DX, stride(DX, 1), DY, stride(DY, 1)) -end - ## nrm2 """ @@ -453,7 +411,10 @@ for (fname, elty, ret_type) in ((:dnrm2_,:Float64,:Float64), end end end -nrm2(x::Union{AbstractVector,DenseArray}) = nrm2(length(x), x, stride1(x)) +# openblas returns 0 for negative stride +nrm2(x::Union{AbstractArray}) = let (p, st) = ptrst1(x) + GC.@preserve x nrm2(length(x), p, abs(st)) +end ## asum @@ -490,8 +451,9 @@ for (fname, elty, ret_type) in ((:dasum_,:Float64,:Float64), end end end -asum(x::Union{AbstractVector,DenseArray}) = asum(length(x), x, stride1(x)) - +asum(x::Union{AbstractArray}) = let (p, st) = ptrst1(x) + GC.@preserve x asum(length(x), p, abs(st)) +end ## axpy """ @@ -538,7 +500,8 @@ function axpy!(alpha::Number, x::Union{DenseArray{T},StridedVector{T}}, y::Union if length(x) != length(y) throw(DimensionMismatch("x has length $(length(x)), but y has length $(length(y))")) end - return axpy!(length(x), convert(T,alpha), x, stride(x, 1), y, stride(y, 1)) + GC.@preserve x y axpy!(length(x), convert(T,alpha), ptrst1(x)..., ptrst1(y)...) + y end function axpy!(alpha::Number, x::Array{T}, rx::Union{UnitRange{Ti},AbstractRange{Ti}}, @@ -555,9 +518,9 @@ function axpy!(alpha::Number, x::Array{T}, rx::Union{UnitRange{Ti},AbstractRange GC.@preserve x y axpy!( length(rx), convert(T, alpha), - pointer(x) + (first(rx) - 1)*sizeof(T), + pointer(x, minimum(rx)), step(rx), - pointer(y) + (first(ry) - 1)*sizeof(T), + pointer(y, minimum(ry)), step(ry)) return y @@ -609,7 +572,8 @@ function axpby!(alpha::Number, x::Union{DenseArray{T},AbstractVector{T}}, beta:: if length(x) != length(y) throw(DimensionMismatch("x has length $(length(x)), but y has length $(length(y))")) end - return axpby!(length(x), convert(T, alpha), x, stride(x, 1), convert(T, beta), y, stride(y, 1)) + GC.@preserve x y axpby!(length(x), convert(T, alpha), ptrst1(x)..., convert(T, beta), ptrst1(y)...) + y end ## iamax @@ -666,10 +630,7 @@ for (fname, elty) in ((:dgemv_,:Float64), chkstride1(A) lda = stride(A,2) lda >= max(1, size(A,1)) || error("`stride(A,2)` must be at least `max(1, size(A,1))`") - sX = stride(X,1) - pX = pointer(X, sX > 0 ? firstindex(X) : lastindex(X)) - sY = stride(Y,1) - pY = pointer(Y, sY > 0 ? firstindex(Y) : lastindex(Y)) + (pX, sX), (pY, sY) = ptrst1(X), ptrst1(Y) GC.@preserve X Y ccall((@blasfunc($fname), libblastrampoline), Cvoid, (Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt}, Ref{$elty}, Ptr{$elty}, Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt}, @@ -750,14 +711,15 @@ for (fname, elty) in ((:dgbmv_,:Float64), y::AbstractVector{$elty}) require_one_based_indexing(A, x, y) chkstride1(A) - ccall((@blasfunc($fname), libblastrampoline), Cvoid, + (px, stx), (py, sty) = ptrst1(x), ptrst1(y) + GC.@preserve x y ccall((@blasfunc($fname), libblastrampoline), Cvoid, (Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt}, Ref{BlasInt}, Ref{BlasInt}, Ref{$elty}, Ptr{$elty}, Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt}, Ref{$elty}, Ptr{$elty}, Ref{BlasInt}, Clong), trans, m, size(A,2), kl, ku, alpha, A, max(1,stride(A,2)), - x, stride(x,1), beta, y, stride(y,1), 1) + px, stx, beta, py, sty, 1) y end function gbmv(trans::AbstractChar, m::Integer, kl::Integer, ku::Integer, alpha::($elty), A::AbstractMatrix{$elty}, x::AbstractVector{$elty}) @@ -810,13 +772,14 @@ for (fname, elty, lib) in ((:dsymv_,:Float64,libblastrampoline), throw(DimensionMismatch("A has size $(size(A)), and y has length $(length(y))")) end chkstride1(A) - ccall((@blasfunc($fname), $lib), Cvoid, + (px, stx), (py, sty) = ptrst1(x), ptrst1(y) + GC.@preserve x y ccall((@blasfunc($fname), $lib), Cvoid, (Ref{UInt8}, Ref{BlasInt}, Ref{$elty}, Ptr{$elty}, Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt}, Ref{$elty}, Ptr{$elty}, Ref{BlasInt}, Clong), uplo, n, alpha, A, - max(1,stride(A,2)), x, stride(x,1), beta, - y, stride(y,1), 1) + max(1,stride(A,2)), px, stx, beta, + py, sty, 1) y end function symv(uplo::AbstractChar, alpha::($elty), A::AbstractMatrix{$elty}, x::AbstractVector{$elty}) @@ -872,15 +835,14 @@ for (fname, elty) in ((:zhemv_,:ComplexF64), end chkstride1(A) lda = max(1, stride(A, 2)) - incx = stride(x, 1) - incy = stride(y, 1) - ccall((@blasfunc($fname), libblastrampoline), Cvoid, + (px, stx), (py, sty) = ptrst1(x), ptrst1(y) + GC.@preserve x y ccall((@blasfunc($fname), libblastrampoline), Cvoid, (Ref{UInt8}, Ref{BlasInt}, Ref{$elty}, Ptr{$elty}, Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt}, Ref{$elty}, Ptr{$elty}, Ref{BlasInt}, Clong), uplo, n, α, A, - lda, x, incx, β, - y, incy, 1) + lda, px, stx, β, + py, sty, 1) y end function hemv(uplo::AbstractChar, α::($elty), A::AbstractMatrix{$elty}, x::AbstractVector{$elty}) @@ -968,7 +930,8 @@ function hpmv!(uplo::AbstractChar, if 2*length(AP) < N*(N + 1) throw(DimensionMismatch("Packed Hermitian matrix A has size smaller than length(x) = $(N).")) end - return hpmv!(uplo, N, convert(T, α), AP, x, stride(x, 1), convert(T, β), y, stride(y, 1)) + GC.@preserve x y hpmv!(uplo, N, convert(T, α), AP, ptrst1(x)..., convert(T, β), ptrst1(y)...) + y end """ @@ -1009,13 +972,14 @@ for (fname, elty) in ((:dsbmv_,:Float64), function sbmv!(uplo::AbstractChar, k::Integer, alpha::($elty), A::AbstractMatrix{$elty}, x::AbstractVector{$elty}, beta::($elty), y::AbstractVector{$elty}) require_one_based_indexing(A, x, y) chkstride1(A) - ccall((@blasfunc($fname), libblastrampoline), Cvoid, + (px, stx), (py, sty) = ptrst1(x), ptrst1(y) + GC.@preserve x y ccall((@blasfunc($fname), libblastrampoline), Cvoid, (Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt}, Ref{$elty}, Ptr{$elty}, Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt}, Ref{$elty}, Ptr{$elty}, Ref{BlasInt}, Clong), uplo, size(A,2), k, alpha, - A, max(1,stride(A,2)), x, stride(x,1), - beta, y, stride(y,1), 1) + A, max(1,stride(A,2)), px, stx, + beta, py, sty, 1) y end function sbmv(uplo::AbstractChar, k::Integer, alpha::($elty), A::AbstractMatrix{$elty}, x::AbstractVector{$elty}) @@ -1118,7 +1082,8 @@ function spmv!(uplo::AbstractChar, if 2*length(AP) < N*(N + 1) throw(DimensionMismatch("Packed symmetric matrix A has size smaller than length(x) = $(N).")) end - return spmv!(uplo, N, convert(T, α), AP, x, stride(x, 1), convert(T, β), y, stride(y, 1)) + GC.@preserve x y spmv!(uplo, N, convert(T, α), AP, ptrst1(x)..., convert(T, β), ptrst1(y)...) + y end """ @@ -1159,13 +1124,14 @@ for (fname, elty) in ((:zhbmv_,:ComplexF64), function hbmv!(uplo::AbstractChar, k::Integer, alpha::($elty), A::AbstractMatrix{$elty}, x::AbstractVector{$elty}, beta::($elty), y::AbstractVector{$elty}) require_one_based_indexing(A, x, y) chkstride1(A) - ccall((@blasfunc($fname), libblastrampoline), Cvoid, + (px, stx), (py, sty) = ptrst1(x), ptrst1(y) + GC.@preserve x y ccall((@blasfunc($fname), libblastrampoline), Cvoid, (Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt}, Ref{$elty}, Ptr{$elty}, Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt}, Ref{$elty}, Ptr{$elty}, Ref{BlasInt}, Clong), uplo, size(A,2), k, alpha, - A, max(1,stride(A,2)), x, stride(x,1), - beta, y, stride(y,1), 1) + A, max(1,stride(A,2)), px, stx, + beta, py, sty, 1) y end function hbmv(uplo::AbstractChar, k::Integer, alpha::($elty), A::AbstractMatrix{$elty}, x::AbstractVector{$elty}) @@ -1219,12 +1185,13 @@ for (fname, elty) in ((:dtrmv_,:Float64), throw(DimensionMismatch("A has size ($n,$n), x has length $(length(x))")) end chkstride1(A) - ccall((@blasfunc($fname), libblastrampoline), Cvoid, + px, stx = ptrst1(x) + GC.@preserve x ccall((@blasfunc($fname), libblastrampoline), Cvoid, (Ref{UInt8}, Ref{UInt8}, Ref{UInt8}, Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt}, Clong, Clong, Clong), uplo, trans, diag, n, - A, max(1,stride(A,2)), x, max(1,stride(x, 1)), 1, 1, 1) + A, max(1,stride(A,2)), px, stx, 1, 1, 1) x end function trmv(uplo::AbstractChar, trans::AbstractChar, diag::AbstractChar, A::AbstractMatrix{$elty}, x::AbstractVector{$elty}) @@ -1274,12 +1241,13 @@ for (fname, elty) in ((:dtrsv_,:Float64), throw(DimensionMismatch("size of A is $n != length(x) = $(length(x))")) end chkstride1(A) - ccall((@blasfunc($fname), libblastrampoline), Cvoid, + px, stx = ptrst1(x) + GC.@preserve x ccall((@blasfunc($fname), libblastrampoline), Cvoid, (Ref{UInt8}, Ref{UInt8}, Ref{UInt8}, Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt}, Clong, Clong, Clong), uplo, trans, diag, n, - A, max(1,stride(A,2)), x, stride(x, 1), 1, 1, 1) + A, max(1,stride(A,2)), px, stx, 1, 1, 1) x end function trsv(uplo::AbstractChar, trans::AbstractChar, diag::AbstractChar, A::AbstractMatrix{$elty}, x::AbstractVector{$elty}) @@ -1993,9 +1961,9 @@ function copyto!(dest::Array{T}, rdest::Union{UnitRange{Ti},AbstractRange{Ti}}, end GC.@preserve src dest BLAS.blascopy!( length(rsrc), - pointer(src) + (first(rsrc) - 1) * sizeof(T), + pointer(src, minimum(rsrc)), step(rsrc), - pointer(dest) + (first(rdest) - 1) * sizeof(T), + pointer(dest, minimum(rdest)), step(rdest)) return dest diff --git a/stdlib/LinearAlgebra/test/blas.jl b/stdlib/LinearAlgebra/test/blas.jl index df29c171b2060..0544b1282c8cc 100644 --- a/stdlib/LinearAlgebra/test/blas.jl +++ b/stdlib/LinearAlgebra/test/blas.jl @@ -70,77 +70,54 @@ Random.seed!(100) end end @testset "iamax" begin - if elty <: Real - x = convert(Vector{elty}, randn(n)) - @test BLAS.iamax(x) == argmax(abs.(x)) - else - z = convert(Vector{elty}, complex.(randn(n),randn(n))) - @test BLAS.iamax(z) == argmax(map(x -> abs(real(x)) + abs(imag(x)), z)) - end + x = convert(Vector{elty}, randn(elty, n)) + @test BLAS.iamax(x) == findmax(x -> abs(real(x)) + abs(imag(x)), x)[2] end @testset "rot!" begin - if elty <: Real - x = convert(Vector{elty}, randn(n)) - y = convert(Vector{elty}, randn(n)) - c = rand(elty) - s = rand(elty) + x = convert(Vector{elty}, randn(elty, n)) + y = convert(Vector{elty}, randn(elty, n)) + c = rand(real(elty)) + for sty in unique!([real(elty), elty]) + s = rand(sty) x2 = copy(x) y2 = copy(y) BLAS.rot!(n, x, 1, y, 1, c, s) @test x ≈ c*x2 + s*y2 - @test y ≈ -s*x2 + c*y2 - else - x = convert(Vector{elty}, complex.(randn(n),rand(n))) - y = convert(Vector{elty}, complex.(randn(n),rand(n))) - cty = (elty == ComplexF32) ? Float32 : Float64 - c = rand(cty) - for sty in [cty, elty] - s = rand(sty) - x2 = copy(x) - y2 = copy(y) - BLAS.rot!(n, x, 1, y, 1, c, s) - @test x ≈ c*x2 + s*y2 - @test y ≈ -conj(s)*x2 + c*y2 - end + @test y ≈ -conj(s)*x2 + c*y2 end end @testset "axp(b)y" begin - if elty <: Real - x1 = convert(Vector{elty}, randn(n)) - x2 = convert(Vector{elty}, randn(n)) - α = rand(elty) - β = rand(elty) - @test BLAS.axpy!(α,copy(x1),copy(x2)) ≈ α*x1 + x2 - @test BLAS.axpby!(α,copy(x1),β,copy(x2)) ≈ α*x1 + β*x2 - @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 - else - z1 = convert(Vector{elty}, complex.(randn(n), randn(n))) - z2 = convert(Vector{elty}, complex.(randn(n), randn(n))) - α = rand(elty) - @test BLAS.axpy!(α, copy(z1), copy(z2)) ≈ z2 + α * z1 - @test_throws DimensionMismatch BLAS.axpy!(α, copy(z1), rand(elty, n + 1)) - @test_throws DimensionMismatch BLAS.axpy!(α, copy(z1), 1:div(n, 2), copy(z2), 1:(div(n, 2) + 1)) - @test_throws ArgumentError BLAS.axpy!(α, copy(z1), 0:div(n,2), copy(z2), 1:(div(n, 2) + 1)) - @test_throws ArgumentError BLAS.axpy!(α, copy(z1), 1:div(n,2), copy(z2), 0:(div(n, 2) - 1)) - @test BLAS.axpy!(α,copy(z1),1:n,copy(z2),1:n) ≈ z2 + α*z1 - end + x1 = convert(Vector{elty}, randn(elty, n)) + x2 = convert(Vector{elty}, 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 + @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) - if elty <: Real - @test BLAS.asum(b) ≈ sum(abs.(b)) - @test BLAS.iamax(b) ≈ argmax(abs.(b)) - else - @test BLAS.asum(b) ≈ sum(abs.(real(b))) + sum(abs.(imag(b))) - @test BLAS.iamax(b) == argmax(map(x -> abs(real(x)) + abs(imag(x)), b)) - end + @test BLAS.asum(b) ≈ sum(x -> abs(real(x)) + abs(imag(x)), b) + @test BLAS.iamax(b) == findmax(x -> abs(real(x)) + abs(imag(x)), b)[2] + # negative stride test + c = view(a,n:-2:2) + @test BLAS.nrm2(c) ≈ norm(c) + @test BLAS.asum(c) ≈ sum(x -> abs(real(x)) + abs(imag(x)), c) + @test BLAS.iamax(c) == 0 end # scal α = rand(elty) @@ -149,8 +126,8 @@ Random.seed!(100) @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)) - @test A\x ≈ BLAS.trsv('U','N','N',A,x) + @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)) end end @@ -185,11 +162,14 @@ Random.seed!(100) @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)) + BLAS.copyto!(x2, 1:n, x1, n:-1:1) + @test x2 == 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 "symmetric/Hermitian multiplication" begin x = rand(elty,n) A = rand(elty,n,n) @@ -197,11 +177,15 @@ Random.seed!(100) 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) 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) end @@ -253,6 +237,8 @@ Random.seed!(100) y_result_blas_lower = copy(y) BLAS.hpmv!('L', α, AP, x, β, y_result_blas_lower) @test y_result_julia_lower ≈ y_result_blas_lower + @test BLAS.hpmv!('L', α, AP, view(x,n:-1:1), β, copy(y)) ≈ + BLAS.hpmv!('L', α, AP, x[n:-1:1], β, copy(y)) y_result_julia_upper = α*AU*x + β*y @@ -267,6 +253,8 @@ Random.seed!(100) y_result_blas_upper = copy(y) BLAS.hpmv!('U', α, AP, x, β, y_result_blas_upper) @test y_result_julia_upper ≈ y_result_blas_upper + @test BLAS.hpmv!('U', α, AP, view(x,n:-1:1), β, copy(y)) ≈ + BLAS.hpmv!('U', α, AP, x[n:-1:1], β, copy(y)) end end @@ -296,7 +284,8 @@ Random.seed!(100) y_result_blas_lower = copy(y) BLAS.spmv!('L', α, AP, x, β, y_result_blas_lower) @test y_result_julia_lower ≈ y_result_blas_lower - + @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 @@ -311,6 +300,8 @@ Random.seed!(100) y_result_blas_upper = copy(y) BLAS.spmv!('U', α, AP, x, β, y_result_blas_upper) @test y_result_julia_upper ≈ y_result_blas_upper + @test BLAS.spmv!('U', α, AP, view(x,n:-1:1), β, copy(y)) ≈ + BLAS.spmv!('U', α, AP, x[n:-1:1], β, copy(y)) end end @@ -330,6 +321,8 @@ Random.seed!(100) 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 end #will work for SymTridiagonal only! @testset "sbmv" begin @@ -341,6 +334,8 @@ Random.seed!(100) 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 else dv = real(rand(elty,n)) ev = rand(elty,n-1) @@ -349,6 +344,8 @@ Random.seed!(100) 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 end end end @@ -552,8 +549,8 @@ end @test BLAS.iamax(x) == 2 M = fill(elty(1.0), 3, 3) - BLAS.scal!(elty(2), view(M,:,2)) - BLAS.scal!(elty(3), view(M,3,:)) + @test BLAS.scal!(elty(2), view(M,:,2)) === view(M,:,2) + @test BLAS.scal!(elty(3), view(M,3,3:-1:1)) === view(M,3,3:-1:1) @test M == elty[1. 2. 1.; 1. 2. 1.; 3. 6. 3.] # Level 2 A = WrappedArray(elty[1 2; 3 4])