diff --git a/base/abstractarray.jl b/base/abstractarray.jl index b119b1d2e356a..e3d0cff63d831 100644 --- a/base/abstractarray.jl +++ b/base/abstractarray.jl @@ -219,46 +219,38 @@ julia> last([1; 2; 3; 4]) last(a) = a[end] """ - stride(A, k::Integer) + strides(A) -Return the distance in memory (in number of elements) between adjacent elements in dimension `k`. +Return a tuple of the memory strides in each dimension. # Examples ```jldoctest julia> A = ones(3,4,5); -julia> stride(A,2) -3 - -julia> stride(A,3) -12 +julia> strides(A) +(1, 3, 12) ``` """ -function stride(a::AbstractArray, i::Integer) - if i > ndims(a) - return length(a) - end - s = 1 - for n = 1:(i-1) - s *= size(a, n) - end - return s -end +function strides end """ - strides(A) + stride(A, k::Integer) -Return a tuple of the memory strides in each dimension. +Return the distance in memory (in number of elements) between adjacent elements in dimension `k`. # Examples ```jldoctest julia> A = ones(3,4,5); -julia> strides(A) -(1, 3, 12) +julia> stride(A,2) +3 + +julia> stride(A,3) +12 ``` """ -strides(A::AbstractArray) = size_to_strides(1, size(A)...) +stride(A::AbstractArray, k::Integer) = strides(A)[k] + @inline size_to_strides(s, d, sz...) = (s, size_to_strides(s * d, sz...)...) size_to_strides(s, d) = (s,) size_to_strides(s) = () diff --git a/base/deprecated.jl b/base/deprecated.jl index 75995dd1d585b..729a8bd45dc25 100644 --- a/base/deprecated.jl +++ b/base/deprecated.jl @@ -3813,6 +3813,21 @@ end @deprecate getq(F::Factorization) F.Q end +# Issues #17812 Remove default stride implementation +function strides(a::AbstractArray) + depwarn(""" + `strides(a::AbstractArray)` is deprecated for general arrays. + Specialize `strides` for custom array types that have the appropriate representation in memory. + Warning: inappropriately implementing this method for an array type that does not use strided + storage may lead to incorrect results or segfaults. + """, :strides) + size_to_strides(1, size(a)...) +end + +@deprecate substrides(s, parent, dim, I::Tuple) substrides(parent, strides(parent), I) + +# END 0.7 deprecations + @deprecate lexcmp(x::AbstractArray, y::AbstractArray) cmp(x, y) @deprecate lexcmp(x::Real, y::Real) cmp(isless, x, y) @deprecate lexcmp(x::Complex, y::Complex) cmp((real(x),imag(x)), (real(y),imag(y))) @@ -3871,6 +3886,7 @@ end @deprecate findin(a, b) find(occursin(b), a) + # END 0.7 deprecations # BEGIN 1.0 deprecations diff --git a/base/linalg/blas.jl b/base/linalg/blas.jl index a31f084c8078e..f65e5f826f9a2 100644 --- a/base/linalg/blas.jl +++ b/base/linalg/blas.jl @@ -63,7 +63,7 @@ export const libblas = Base.libblas_name const liblapack = Base.liblapack_name -import ..LinAlg: BlasReal, BlasComplex, BlasFloat, BlasInt, DimensionMismatch, checksquare, axpy! +import ..LinAlg: BlasReal, BlasComplex, BlasFloat, BlasInt, DimensionMismatch, checksquare, stride1, chkstride1, axpy! # utility routines function vendor() @@ -177,7 +177,7 @@ for (fname, elty) in ((:dcopy_,:Float64), (:ccopy_,:ComplexF32)) @eval begin # SUBROUTINE DCOPY(N,DX,INCX,DY,INCY) - function blascopy!(n::Integer, DX::Union{Ptr{$elty},StridedArray{$elty}}, incx::Integer, DY::Union{Ptr{$elty},StridedArray{$elty}}, incy::Integer) + function blascopy!(n::Integer, DX::Union{Ptr{$elty},AbstractArray{$elty}}, incx::Integer, DY::Union{Ptr{$elty},AbstractArray{$elty}}, incy::Integer) ccall((@blasfunc($fname), libblas), Cvoid, (Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt}), n, DX, incx, DY, incy) @@ -208,7 +208,7 @@ for (fname, elty) in ((:dscal_,:Float64), (:cscal_,:ComplexF32)) @eval begin # SUBROUTINE DSCAL(N,DA,DX,INCX) - function scal!(n::Integer, DA::$elty, DX::Union{Ptr{$elty},DenseArray{$elty}}, incx::Integer) + function scal!(n::Integer, DA::$elty, DX::Union{Ptr{$elty},AbstractArray{$elty}}, incx::Integer) ccall((@blasfunc($fname), libblas), Cvoid, (Ref{BlasInt}, Ref{$elty}, Ptr{$elty}, Ref{BlasInt}), n, DA, DX, incx) @@ -272,7 +272,7 @@ for (fname, elty) in ((:ddot_,:Float64), # * .. # * .. Array Arguments .. # DOUBLE PRECISION DX(*),DY(*) - function dot(n::Integer, DX::Union{Ptr{$elty},DenseArray{$elty}}, incx::Integer, DY::Union{Ptr{$elty},DenseArray{$elty}}, incy::Integer) + function dot(n::Integer, DX::Union{Ptr{$elty},AbstractArray{$elty}}, incx::Integer, DY::Union{Ptr{$elty},AbstractArray{$elty}}, incy::Integer) ccall((@blasfunc($fname), libblas), $elty, (Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt}), n, DX, incx, DY, incy) @@ -288,7 +288,7 @@ for (fname, elty) in ((:cblas_zdotc_sub,:ComplexF64), # * .. # * .. Array Arguments .. # DOUBLE PRECISION DX(*),DY(*) - function dotc(n::Integer, DX::Union{Ptr{$elty},DenseArray{$elty}}, incx::Integer, DY::Union{Ptr{$elty},DenseArray{$elty}}, incy::Integer) + function dotc(n::Integer, DX::Union{Ptr{$elty},AbstractArray{$elty}}, incx::Integer, DY::Union{Ptr{$elty},AbstractArray{$elty}}, incy::Integer) result = Ref{$elty}() ccall((@blasfunc($fname), libblas), Cvoid, (BlasInt, Ptr{$elty}, BlasInt, Ptr{$elty}, BlasInt, Ptr{$elty}), @@ -306,7 +306,7 @@ for (fname, elty) in ((:cblas_zdotu_sub,:ComplexF64), # * .. # * .. Array Arguments .. # DOUBLE PRECISION DX(*),DY(*) - function dotu(n::Integer, DX::Union{Ptr{$elty},DenseArray{$elty}}, incx::Integer, DY::Union{Ptr{$elty},DenseArray{$elty}}, incy::Integer) + function dotu(n::Integer, DX::Union{Ptr{$elty},AbstractArray{$elty}}, incx::Integer, DY::Union{Ptr{$elty},AbstractArray{$elty}}, incy::Integer) result = Ref{$elty}() ccall((@blasfunc($fname), libblas), Cvoid, (BlasInt, Ptr{$elty}, BlasInt, Ptr{$elty}, BlasInt, Ptr{$elty}), @@ -315,21 +315,22 @@ for (fname, elty) in ((:cblas_zdotu_sub,:ComplexF64), end end end -function dot(DX::Union{DenseArray{T},StridedVector{T}}, DY::Union{DenseArray{T},StridedVector{T}}) where T<:BlasReal + +function dot(DX::Union{DenseArray{T},AbstractVector{T}}, DY::Union{DenseArray{T},AbstractVector{T}}) where T<:BlasReal n = length(DX) if n != length(DY) throw(DimensionMismatch("dot product arguments have lengths $(length(DX)) and $(length(DY))")) end Base.@gc_preserve DX DY dot(n, pointer(DX), stride(DX, 1), pointer(DY), stride(DY, 1)) end -function dotc(DX::Union{DenseArray{T},StridedVector{T}}, DY::Union{DenseArray{T},StridedVector{T}}) where T<:BlasComplex +function dotc(DX::Union{DenseArray{T},AbstractVector{T}}, DY::Union{DenseArray{T},AbstractVector{T}}) where T<:BlasComplex n = length(DX) if n != length(DY) throw(DimensionMismatch("dot product arguments have lengths $(length(DX)) and $(length(DY))")) end Base.@gc_preserve DX DY dotc(n, pointer(DX), stride(DX, 1), pointer(DY), stride(DY, 1)) end -function dotu(DX::Union{DenseArray{T},StridedVector{T}}, DY::Union{DenseArray{T},StridedVector{T}}) where T<:BlasComplex +function dotu(DX::Union{DenseArray{T},AbstractVector{T}}, DY::Union{DenseArray{T},AbstractVector{T}}) where T<:BlasComplex n = length(DX) if n != length(DY) throw(DimensionMismatch("dot product arguments have lengths $(length(DX)) and $(length(DY))")) @@ -339,9 +340,6 @@ end ## nrm2 -stride1(x) = stride(x,1) -stride1(x::Array) = 1 - """ nrm2(n, X, incx) @@ -364,14 +362,14 @@ for (fname, elty, ret_type) in ((:dnrm2_,:Float64,:Float64), (:scnrm2_,:ComplexF32,:Float32)) @eval begin # SUBROUTINE DNRM2(N,X,INCX) - function nrm2(n::Integer, X::Union{Ptr{$elty},DenseArray{$elty}}, incx::Integer) + function nrm2(n::Integer, X::Union{Ptr{$elty},AbstractArray{$elty}}, incx::Integer) ccall((@blasfunc($fname), libblas), $ret_type, (Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt}), n, X, incx) end end end -nrm2(x::Union{StridedVector,Array}) = Base.@gc_preserve x nrm2(length(x), pointer(x), stride1(x)) +nrm2(x::Union{AbstractVector,DenseArray}) = Base.@gc_preserve x nrm2(length(x), pointer(x), stride1(x)) ## asum @@ -397,14 +395,14 @@ for (fname, elty, ret_type) in ((:dasum_,:Float64,:Float64), (:scasum_,:ComplexF32,:Float32)) @eval begin # SUBROUTINE ASUM(N, X, INCX) - function asum(n::Integer, X::Union{Ptr{$elty},DenseArray{$elty}}, incx::Integer) + function asum(n::Integer, X::Union{Ptr{$elty},AbstractArray{$elty}}, incx::Integer) ccall((@blasfunc($fname), libblas), $ret_type, (Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt}), n, X, incx) end end end -asum(x::Union{StridedVector,Array}) = Base.@gc_preserve x asum(length(x), pointer(x), stride1(x)) +asum(x::Union{AbstractVector,DenseArray}) = Base.@gc_preserve x asum(length(x), pointer(x), stride1(x)) ## axpy @@ -440,7 +438,7 @@ for (fname, elty) in ((:daxpy_,:Float64), # INTEGER INCX,INCY,N #* .. Array Arguments .. # DOUBLE PRECISION DX(*),DY(*) - function axpy!(n::Integer, alpha::($elty), dx::Union{Ptr{$elty}, DenseArray{$elty}}, incx::Integer, dy::Union{Ptr{$elty}, DenseArray{$elty}}, incy::Integer) + function axpy!(n::Integer, alpha::($elty), dx::Union{Ptr{$elty}, AbstractArray{$elty}}, incx::Integer, dy::Union{Ptr{$elty}, AbstractArray{$elty}}, incy::Integer) ccall((@blasfunc($fname), libblas), Cvoid, (Ref{BlasInt}, Ref{$elty}, Ptr{$elty}, Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt}), n, alpha, dx, incx, dy, incy) @@ -502,8 +500,8 @@ for (fname, elty) in ((:daxpby_,:Float64), (:saxpby_,:Float32), #* .. Array Arguments .. # DOUBLE PRECISION DX(*),DY(*) function axpby!(n::Integer, alpha::($elty), dx::Union{Ptr{$elty}, - DenseArray{$elty}}, incx::Integer, beta::($elty), - dy::Union{Ptr{$elty}, DenseArray{$elty}}, incy::Integer) + AbstractArray{$elty}}, incx::Integer, beta::($elty), + dy::Union{Ptr{$elty}, AbstractArray{$elty}}, incy::Integer) ccall((@blasfunc($fname), libblas), Cvoid, (Ref{BlasInt}, Ref{$elty}, Ptr{$elty}, Ref{BlasInt}, Ref{$elty}, Ptr{$elty}, Ref{BlasInt}), n, alpha, dx, incx, beta, dy, incy) @@ -512,7 +510,7 @@ for (fname, elty) in ((:daxpby_,:Float64), (:saxpby_,:Float32), end end -function axpby!(alpha::Number, x::Union{DenseArray{T},StridedVector{T}}, beta::Number, y::Union{DenseArray{T},StridedVector{T}}) where T<:BlasFloat +function axpby!(alpha::Number, x::Union{DenseArray{T},AbstractVector{T}}, beta::Number, y::Union{DenseArray{T},AbstractVector{T}}) where T<:BlasFloat if length(x) != length(y) throw(DimensionMismatch("x has length $(length(x)), but y has length $(length(y))")) end @@ -526,14 +524,14 @@ for (fname, elty) in ((:idamax_,:Float64), (:izamax_,:ComplexF64), (:icamax_,:ComplexF32)) @eval begin - function iamax(n::Integer, dx::Union{Ptr{$elty}, DenseArray{$elty}}, incx::Integer) + function iamax(n::Integer, dx::Union{Ptr{$elty}, AbstractArray{$elty}}, incx::Integer) ccall((@blasfunc($fname), libblas),BlasInt, (Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt}), n, dx, incx) end end end -iamax(dx::Union{StridedVector,Array}) = Base.@gc_preserve dx iamax(length(dx), pointer(dx), stride1(dx)) +iamax(dx::Union{AbstractVector,DenseArray}) = Base.@gc_preserve dx iamax(length(dx), pointer(dx), stride1(dx)) # Level 2 ## mv @@ -550,7 +548,7 @@ for (fname, elty) in ((:dgemv_,:Float64), # CHARACTER TRANS #* .. Array Arguments .. # DOUBLE PRECISION A(LDA,*),X(*),Y(*) - function gemv!(trans::Char, alpha::($elty), A::StridedVecOrMat{$elty}, X::StridedVector{$elty}, beta::($elty), Y::StridedVector{$elty}) + function gemv!(trans::Char, alpha::($elty), A::AbstractVecOrMat{$elty}, X::AbstractVector{$elty}, beta::($elty), Y::AbstractVector{$elty}) m,n = size(A,1),size(A,2) if trans == 'N' && (length(X) != n || length(Y) != m) throw(DimensionMismatch("A has dimensions $(size(A)), X has length $(length(X)) and Y has length $(length(Y))")) @@ -559,6 +557,7 @@ for (fname, elty) in ((:dgemv_,:Float64), elseif trans == 'T' && (length(X) != m || length(Y) != n) throw(DimensionMismatch("the transpose of A has dimensions $n, $m, X has length $(length(X)) and Y has length $(length(Y))")) end + chkstride1(A) ccall((@blasfunc($fname), libblas), Cvoid, (Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt}, Ref{$elty}, Ptr{$elty}, Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt}, @@ -568,10 +567,10 @@ for (fname, elty) in ((:dgemv_,:Float64), beta, Y, stride(Y,1)) Y end - function gemv(trans::Char, alpha::($elty), A::StridedMatrix{$elty}, X::StridedVector{$elty}) + function gemv(trans::Char, alpha::($elty), A::AbstractMatrix{$elty}, X::AbstractVector{$elty}) gemv!(trans, alpha, A, X, zero($elty), similar(X, $elty, size(A, (trans == 'N' ? 1 : 2)))) end - function gemv(trans::Char, A::StridedMatrix{$elty}, X::StridedVector{$elty}) + function gemv(trans::Char, A::AbstractMatrix{$elty}, X::AbstractVector{$elty}) gemv!(trans, one($elty), A, X, zero($elty), similar(X, $elty, size(A, (trans == 'N' ? 1 : 2)))) end end @@ -633,7 +632,8 @@ for (fname, elty) in ((:dgbmv_,:Float64), # CHARACTER TRANS # * .. Array Arguments .. # DOUBLE PRECISION A(LDA,*),X(*),Y(*) - function gbmv!(trans::Char, m::Integer, kl::Integer, ku::Integer, alpha::($elty), A::StridedMatrix{$elty}, x::StridedVector{$elty}, beta::($elty), y::StridedVector{$elty}) + function gbmv!(trans::Char, m::Integer, kl::Integer, ku::Integer, alpha::($elty), A::AbstractMatrix{$elty}, x::AbstractVector{$elty}, beta::($elty), y::AbstractVector{$elty}) + chkstride1(A) ccall((@blasfunc($fname), libblas), Cvoid, (Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt}, Ref{BlasInt}, Ref{BlasInt}, Ref{$elty}, Ptr{$elty}, Ref{BlasInt}, @@ -644,12 +644,12 @@ for (fname, elty) in ((:dgbmv_,:Float64), x, stride(x,1), beta, y, stride(y,1)) y end - function gbmv(trans::Char, m::Integer, kl::Integer, ku::Integer, alpha::($elty), A::StridedMatrix{$elty}, x::StridedVector{$elty}) + function gbmv(trans::Char, m::Integer, kl::Integer, ku::Integer, alpha::($elty), A::AbstractMatrix{$elty}, x::AbstractVector{$elty}) n = size(A,2) leny = trans == 'N' ? m : n gbmv!(trans, m, kl, ku, alpha, A, x, zero($elty), similar(x, $elty, leny)) end - function gbmv(trans::Char, m::Integer, kl::Integer, ku::Integer, A::StridedMatrix{$elty}, x::StridedVector{$elty}) + function gbmv(trans::Char, m::Integer, kl::Integer, ku::Integer, A::AbstractMatrix{$elty}, x::AbstractVector{$elty}) gbmv(trans, m, kl, ku, one($elty), A, x) end end @@ -679,7 +679,7 @@ for (fname, elty, lib) in ((:dsymv_,:Float64,libblas), # CHARACTER UPLO # .. Array Arguments .. # DOUBLE PRECISION A(LDA,*),X(*),Y(*) - function symv!(uplo::Char, alpha::($elty), A::StridedMatrix{$elty}, x::StridedVector{$elty},beta::($elty), y::StridedVector{$elty}) + function symv!(uplo::Char, alpha::($elty), A::AbstractMatrix{$elty}, x::AbstractVector{$elty}, beta::($elty), y::AbstractVector{$elty}) m, n = size(A) if m != n throw(DimensionMismatch("matrix A is $m by $n but must be square")) @@ -690,6 +690,7 @@ for (fname, elty, lib) in ((:dsymv_,:Float64,libblas), if m != length(y) throw(DimensionMismatch("A has size $(size(A)), and y has length $(length(y))")) end + chkstride1(A) ccall((@blasfunc($fname), $lib), Cvoid, (Ref{UInt8}, Ref{BlasInt}, Ref{$elty}, Ptr{$elty}, Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt}, Ref{$elty}, @@ -699,10 +700,10 @@ for (fname, elty, lib) in ((:dsymv_,:Float64,libblas), y, stride(y,1)) y end - function symv(uplo::Char, alpha::($elty), A::StridedMatrix{$elty}, x::StridedVector{$elty}) + function symv(uplo::Char, alpha::($elty), A::AbstractMatrix{$elty}, x::AbstractVector{$elty}) symv!(uplo, alpha, A, x, zero($elty), similar(x)) end - function symv(uplo::Char, A::StridedMatrix{$elty}, x::StridedVector{$elty}) + function symv(uplo::Char, A::AbstractMatrix{$elty}, x::AbstractVector{$elty}) symv(uplo, one($elty), A, x) end end @@ -729,7 +730,7 @@ symv(ul, A, x) for (fname, elty) in ((:zhemv_,:ComplexF64), (:chemv_,:ComplexF32)) @eval begin - function hemv!(uplo::Char, α::$elty, A::StridedMatrix{$elty}, x::StridedVector{$elty}, β::$elty, y::StridedVector{$elty}) + function hemv!(uplo::Char, α::$elty, A::AbstractMatrix{$elty}, x::AbstractVector{$elty}, β::$elty, y::AbstractVector{$elty}) m, n = size(A) if m != n throw(DimensionMismatch("matrix A is $m by $n but must be square")) @@ -740,6 +741,7 @@ for (fname, elty) in ((:zhemv_,:ComplexF64), if m != length(y) throw(DimensionMismatch("A has size $(size(A)), and y has length $(length(y))")) end + chkstride1(A) lda = max(1, stride(A, 2)) incx = stride(x, 1) incy = stride(y, 1) @@ -752,10 +754,10 @@ for (fname, elty) in ((:zhemv_,:ComplexF64), y, incy) y end - function hemv(uplo::Char, α::($elty), A::StridedMatrix{$elty}, x::StridedVector{$elty}) + function hemv(uplo::Char, α::($elty), A::AbstractMatrix{$elty}, x::AbstractVector{$elty}) hemv!(uplo, α, A, x, zero($elty), similar(x)) end - function hemv(uplo::Char, A::StridedMatrix{$elty}, x::StridedVector{$elty}) + function hemv(uplo::Char, A::AbstractMatrix{$elty}, x::AbstractVector{$elty}) hemv(uplo, one($elty), A, x) end end @@ -772,7 +774,8 @@ for (fname, elty) in ((:dsbmv_,:Float64), # CHARACTER UPLO # * .. Array Arguments .. # DOUBLE PRECISION A(LDA,*),X(*),Y(*) - function sbmv!(uplo::Char, k::Integer, alpha::($elty), A::StridedMatrix{$elty}, x::StridedVector{$elty}, beta::($elty), y::StridedVector{$elty}) + function sbmv!(uplo::Char, k::Integer, alpha::($elty), A::AbstractMatrix{$elty}, x::AbstractVector{$elty}, beta::($elty), y::AbstractVector{$elty}) + chkstride1(A) ccall((@blasfunc($fname), libblas), Cvoid, (Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt}, Ref{$elty}, Ptr{$elty}, Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt}, @@ -782,11 +785,11 @@ for (fname, elty) in ((:dsbmv_,:Float64), beta, y, stride(y,1)) y end - function sbmv(uplo::Char, k::Integer, alpha::($elty), A::StridedMatrix{$elty}, x::StridedVector{$elty}) + function sbmv(uplo::Char, k::Integer, alpha::($elty), A::AbstractMatrix{$elty}, x::AbstractVector{$elty}) n = size(A,2) sbmv!(uplo, k, alpha, A, x, zero($elty), similar(x, $elty, n)) end - function sbmv(uplo::Char, k::Integer, A::StridedMatrix{$elty}, x::StridedVector{$elty}) + function sbmv(uplo::Char, k::Integer, A::AbstractMatrix{$elty}, x::AbstractVector{$elty}) sbmv(uplo, k, one($elty), A, x) end end @@ -834,7 +837,8 @@ for (fname, elty) in ((:zhbmv_,:ComplexF64), # CHARACTER UPLO # * .. Array Arguments .. # DOUBLE PRECISION A(LDA,*),X(*),Y(*) - function hbmv!(uplo::Char, k::Integer, alpha::($elty), A::StridedMatrix{$elty}, x::StridedVector{$elty}, beta::($elty), y::StridedVector{$elty}) + function hbmv!(uplo::Char, k::Integer, alpha::($elty), A::AbstractMatrix{$elty}, x::AbstractVector{$elty}, beta::($elty), y::AbstractVector{$elty}) + chkstride1(A) ccall((@blasfunc($fname), libblas), Cvoid, (Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt}, Ref{$elty}, Ptr{$elty}, Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt}, @@ -844,11 +848,11 @@ for (fname, elty) in ((:zhbmv_,:ComplexF64), beta, y, stride(y,1)) y end - function hbmv(uplo::Char, k::Integer, alpha::($elty), A::StridedMatrix{$elty}, x::StridedVector{$elty}) + function hbmv(uplo::Char, k::Integer, alpha::($elty), A::AbstractMatrix{$elty}, x::AbstractVector{$elty}) n = size(A,2) hbmv!(uplo, k, alpha, A, x, zero($elty), similar(x, $elty, n)) end - function hbmv(uplo::Char, k::Integer, A::StridedMatrix{$elty}, x::StridedVector{$elty}) + function hbmv(uplo::Char, k::Integer, A::AbstractMatrix{$elty}, x::AbstractVector{$elty}) hbmv(uplo, k, one($elty), A, x) end end @@ -888,11 +892,12 @@ for (fname, elty) in ((:dtrmv_,:Float64), # CHARACTER DIAG,TRANS,UPLO # * .. Array Arguments .. # DOUBLE PRECISION A(LDA,*),X(*) - function trmv!(uplo::Char, trans::Char, diag::Char, A::StridedMatrix{$elty}, x::StridedVector{$elty}) + function trmv!(uplo::Char, trans::Char, diag::Char, A::AbstractMatrix{$elty}, x::AbstractVector{$elty}) n = checksquare(A) if n != length(x) throw(DimensionMismatch("A has size ($n,$n), x has length $(length(x))")) end + chkstride1(A) ccall((@blasfunc($fname), libblas), Cvoid, (Ref{UInt8}, Ref{UInt8}, Ref{UInt8}, Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt}), @@ -900,7 +905,7 @@ for (fname, elty) in ((:dtrmv_,:Float64), A, max(1,stride(A,2)), x, max(1,stride(x, 1))) x end - function trmv(uplo::Char, trans::Char, diag::Char, A::StridedMatrix{$elty}, x::StridedVector{$elty}) + function trmv(uplo::Char, trans::Char, diag::Char, A::AbstractMatrix{$elty}, x::AbstractVector{$elty}) trmv!(uplo, trans, diag, A, copy(x)) end end @@ -940,11 +945,12 @@ for (fname, elty) in ((:dtrsv_,:Float64), # CHARACTER DIAG,TRANS,UPLO # .. Array Arguments .. # DOUBLE PRECISION A(LDA,*),X(*) - function trsv!(uplo::Char, trans::Char, diag::Char, A::StridedMatrix{$elty}, x::StridedVector{$elty}) + function trsv!(uplo::Char, trans::Char, diag::Char, A::AbstractMatrix{$elty}, x::AbstractVector{$elty}) n = checksquare(A) if n != length(x) throw(DimensionMismatch("size of A is $n != length(x) = $(length(x))")) end + chkstride1(A) ccall((@blasfunc($fname), libblas), Cvoid, (Ref{UInt8}, Ref{UInt8}, Ref{UInt8}, Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt}), @@ -952,7 +958,7 @@ for (fname, elty) in ((:dtrsv_,:Float64), A, max(1,stride(A,2)), x, stride(x, 1)) x end - function trsv(uplo::Char, trans::Char, diag::Char, A::StridedMatrix{$elty}, x::StridedVector{$elty}) + function trsv(uplo::Char, trans::Char, diag::Char, A::AbstractMatrix{$elty}, x::AbstractVector{$elty}) trsv!(uplo, trans, diag, A, copy(x)) end end @@ -972,7 +978,7 @@ for (fname, elty) in ((:dger_,:Float64), (:zgerc_,:ComplexF64), (:cgerc_,:ComplexF32)) @eval begin - function ger!(α::$elty, x::StridedVector{$elty}, y::StridedVector{$elty}, A::StridedMatrix{$elty}) + function ger!(α::$elty, x::AbstractVector{$elty}, y::AbstractVector{$elty}, A::AbstractMatrix{$elty}) m, n = size(A) if m != length(x) || n != length(y) throw(DimensionMismatch("A has size ($m,$n), x has length $(length(x)), y has length $(length(y))")) @@ -1004,7 +1010,7 @@ for (fname, elty, lib) in ((:dsyr_,:Float64,libblas), (:zsyr_,:ComplexF64,liblapack), (:csyr_,:ComplexF32,liblapack)) @eval begin - function syr!(uplo::Char, α::$elty, x::StridedVector{$elty}, A::StridedMatrix{$elty}) + function syr!(uplo::Char, α::$elty, x::AbstractVector{$elty}, A::AbstractMatrix{$elty}) n = checksquare(A) if length(x) != n throw(DimensionMismatch("A has size ($n,$n), x has length $(length(x))")) @@ -1033,7 +1039,7 @@ function her! end for (fname, elty, relty) in ((:zher_,:ComplexF64, :Float64), (:cher_,:ComplexF32, :Float32)) @eval begin - function her!(uplo::Char, α::$relty, x::StridedVector{$elty}, A::StridedMatrix{$elty}) + function her!(uplo::Char, α::$relty, x::AbstractVector{$elty}, A::AbstractMatrix{$elty}) n = checksquare(A) if length(x) != n throw(DimensionMismatch("A has size ($n,$n), x has length $(length(x))")) @@ -1072,7 +1078,7 @@ for (gemm, elty) in # CHARACTER TRANSA,TRANSB # * .. Array Arguments .. # DOUBLE PRECISION A(LDA,*),B(LDB,*),C(LDC,*) - function gemm!(transA::Char, transB::Char, alpha::($elty), A::StridedVecOrMat{$elty}, B::StridedVecOrMat{$elty}, beta::($elty), C::StridedVecOrMat{$elty}) + function gemm!(transA::Char, transB::Char, alpha::($elty), A::AbstractVecOrMat{$elty}, B::AbstractVecOrMat{$elty}, beta::($elty), C::AbstractVecOrMat{$elty}) # if any([stride(A,1), stride(B,1), stride(C,1)] .!= 1) # error("gemm!: BLAS module requires contiguous matrix columns") # end # should this be checked on every call? @@ -1083,6 +1089,9 @@ for (gemm, elty) in if ka != kb || m != size(C,1) || n != size(C,2) throw(DimensionMismatch("A has size ($m,$ka), B has size ($kb,$n), C has size $(size(C))")) end + chkstride1(A) + chkstride1(B) + chkstride1(C) ccall((@blasfunc($gemm), libblas), Cvoid, (Ref{UInt8}, Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt}, Ref{BlasInt}, Ref{$elty}, Ptr{$elty}, Ref{BlasInt}, @@ -1094,10 +1103,10 @@ for (gemm, elty) in max(1,stride(C,2))) C end - function gemm(transA::Char, transB::Char, alpha::($elty), A::StridedMatrix{$elty}, B::StridedMatrix{$elty}) + function gemm(transA::Char, transB::Char, alpha::($elty), A::AbstractMatrix{$elty}, B::AbstractMatrix{$elty}) gemm!(transA, transB, alpha, A, B, zero($elty), similar(B, $elty, (size(A, transA == 'N' ? 1 : 2), size(B, transB == 'N' ? 2 : 1)))) end - function gemm(transA::Char, transB::Char, A::StridedMatrix{$elty}, B::StridedMatrix{$elty}) + function gemm(transA::Char, transB::Char, A::AbstractMatrix{$elty}, B::AbstractMatrix{$elty}) gemm(transA, transB, one($elty), A, B) end end @@ -1131,7 +1140,7 @@ for (mfname, elty) in ((:dsymm_,:Float64), # CHARACTER SIDE,UPLO # .. Array Arguments .. # DOUBLE PRECISION A(LDA,*),B(LDB,*),C(LDC,*) - function symm!(side::Char, uplo::Char, alpha::($elty), A::StridedMatrix{$elty}, B::StridedMatrix{$elty}, beta::($elty), C::StridedMatrix{$elty}) + function symm!(side::Char, uplo::Char, alpha::($elty), A::AbstractMatrix{$elty}, B::AbstractMatrix{$elty}, beta::($elty), C::AbstractMatrix{$elty}) m, n = size(C) j = checksquare(A) if j != (side == 'L' ? m : n) @@ -1140,6 +1149,9 @@ for (mfname, elty) in ((:dsymm_,:Float64), if size(B,2) != n throw(DimensionMismatch("B has second dimension $(size(B,2)) but needs to match second dimension of C, $n")) end + chkstride1(A) + chkstride1(B) + chkstride1(C) ccall((@blasfunc($mfname), libblas), Cvoid, (Ref{UInt8}, Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt}, Ref{$elty}, Ptr{$elty}, Ref{BlasInt}, Ptr{$elty}, @@ -1149,10 +1161,10 @@ for (mfname, elty) in ((:dsymm_,:Float64), max(1,stride(B,2)), beta, C, max(1,stride(C,2))) C end - function symm(side::Char, uplo::Char, alpha::($elty), A::StridedMatrix{$elty}, B::StridedMatrix{$elty}) + function symm(side::Char, uplo::Char, alpha::($elty), A::AbstractMatrix{$elty}, B::AbstractMatrix{$elty}) symm!(side, uplo, alpha, A, B, zero($elty), similar(B)) end - function symm(side::Char, uplo::Char, A::StridedMatrix{$elty}, B::StridedMatrix{$elty}) + function symm(side::Char, uplo::Char, A::AbstractMatrix{$elty}, B::AbstractMatrix{$elty}) symm(side, uplo, one($elty), A, B) end end @@ -1196,7 +1208,7 @@ for (mfname, elty) in ((:zhemm_,:ComplexF64), # CHARACTER SIDE,UPLO # .. Array Arguments .. # DOUBLE PRECISION A(LDA,*),B(LDB,*),C(LDC,*) - function hemm!(side::Char, uplo::Char, alpha::($elty), A::StridedMatrix{$elty}, B::StridedMatrix{$elty}, beta::($elty), C::StridedMatrix{$elty}) + function hemm!(side::Char, uplo::Char, alpha::($elty), A::AbstractMatrix{$elty}, B::AbstractMatrix{$elty}, beta::($elty), C::AbstractMatrix{$elty}) m, n = size(C) j = checksquare(A) if j != (side == 'L' ? m : n) @@ -1205,6 +1217,9 @@ for (mfname, elty) in ((:zhemm_,:ComplexF64), if size(B,2) != n throw(DimensionMismatch("B has second dimension $(size(B,2)) but needs to match second dimension of C, $n")) end + chkstride1(A) + chkstride1(B) + chkstride1(C) ccall((@blasfunc($mfname), libblas), Cvoid, (Ref{UInt8}, Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt}, Ref{$elty}, Ptr{$elty}, Ref{BlasInt}, Ptr{$elty}, @@ -1214,10 +1229,10 @@ for (mfname, elty) in ((:zhemm_,:ComplexF64), max(1,stride(B,2)), beta, C, max(1,stride(C,2))) C end - function hemm(side::Char, uplo::Char, alpha::($elty), A::StridedMatrix{$elty}, B::StridedMatrix{$elty}) + function hemm(side::Char, uplo::Char, alpha::($elty), A::AbstractMatrix{$elty}, B::AbstractMatrix{$elty}) hemm!(side, uplo, alpha, A, B, zero($elty), similar(B)) end - function hemm(side::Char, uplo::Char, A::StridedMatrix{$elty}, B::StridedMatrix{$elty}) + function hemm(side::Char, uplo::Char, A::AbstractMatrix{$elty}, B::AbstractMatrix{$elty}) hemm(side, uplo, one($elty), A, B) end end @@ -1257,12 +1272,14 @@ for (fname, elty) in ((:dsyrk_,:Float64), # * .. Array Arguments .. # REAL A(LDA,*),C(LDC,*) function syrk!(uplo::Char, trans::Char, - alpha::($elty), A::StridedVecOrMat{$elty}, - beta::($elty), C::StridedMatrix{$elty}) + alpha::($elty), A::AbstractVecOrMat{$elty}, + beta::($elty), C::AbstractMatrix{$elty}) n = checksquare(C) nn = size(A, trans == 'N' ? 1 : 2) if nn != n throw(DimensionMismatch("C has size ($n,$n), corresponding dimension of A is $nn")) end k = size(A, trans == 'N' ? 2 : 1) + chkstride1(A) + chkstride1(C) ccall((@blasfunc($fname), libblas), Cvoid, (Ref{UInt8}, Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt}, Ref{$elty}, Ptr{$elty}, Ref{BlasInt}, Ref{$elty}, @@ -1274,12 +1291,12 @@ for (fname, elty) in ((:dsyrk_,:Float64), end end end -function syrk(uplo::Char, trans::Char, alpha::Number, A::StridedVecOrMat) +function syrk(uplo::Char, trans::Char, alpha::Number, A::AbstractVecOrMat) T = eltype(A) n = size(A, trans == 'N' ? 1 : 2) syrk!(uplo, trans, convert(T,alpha), A, zero(T), similar(A, T, (n, n))) end -syrk(uplo::Char, trans::Char, A::StridedVecOrMat) = syrk(uplo, trans, one(eltype(A)), A) +syrk(uplo::Char, trans::Char, A::AbstractVecOrMat) = syrk(uplo, trans, one(eltype(A)), A) """ herk!(uplo, trans, alpha, A, beta, C) @@ -1311,13 +1328,15 @@ for (fname, elty, relty) in ((:zherk_, :ComplexF64, :Float64), # * .. # * .. Array Arguments .. # COMPLEX A(LDA,*),C(LDC,*) - function herk!(uplo::Char, trans::Char, α::$relty, A::StridedVecOrMat{$elty}, - β::$relty, C::StridedMatrix{$elty}) + function herk!(uplo::Char, trans::Char, α::$relty, A::AbstractVecOrMat{$elty}, + β::$relty, C::AbstractMatrix{$elty}) n = checksquare(C) nn = size(A, trans == 'N' ? 1 : 2) if nn != n throw(DimensionMismatch("the matrix to update has dimension $n but the implied dimension of the update is $(size(A, trans == 'N' ? 1 : 2))")) end + chkstride1(A) + chkstride1(C) k = size(A, trans == 'N' ? 2 : 1) ccall((@blasfunc($fname), libblas), Cvoid, (Ref{UInt8}, Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt}, @@ -1328,11 +1347,11 @@ for (fname, elty, relty) in ((:zherk_, :ComplexF64, :Float64), C, max(1,stride(C,2))) C end - function herk(uplo::Char, trans::Char, α::$relty, A::StridedVecOrMat{$elty}) + function herk(uplo::Char, trans::Char, α::$relty, A::AbstractVecOrMat{$elty}) n = size(A, trans == 'N' ? 1 : 2) herk!(uplo, trans, α, A, zero($relty), similar(A, (n,n))) end - herk(uplo::Char, trans::Char, A::StridedVecOrMat{$elty}) = herk(uplo, trans, one($relty), A) + herk(uplo::Char, trans::Char, A::AbstractVecOrMat{$elty}) = herk(uplo, trans, one($relty), A) end end @@ -1352,12 +1371,15 @@ for (fname, elty) in ((:dsyr2k_,:Float64), # .. Array Arguments .. # REAL PRECISION A(LDA,*),B(LDB,*),C(LDC,*) function syr2k!(uplo::Char, trans::Char, - alpha::($elty), A::StridedVecOrMat{$elty}, B::StridedVecOrMat{$elty}, - beta::($elty), C::StridedMatrix{$elty}) + alpha::($elty), A::AbstractVecOrMat{$elty}, B::AbstractVecOrMat{$elty}, + beta::($elty), C::AbstractMatrix{$elty}) n = checksquare(C) nn = size(A, trans == 'N' ? 1 : 2) if nn != n throw(DimensionMismatch("C has size ($n,$n), corresponding dimension of A is $nn")) end k = size(A, trans == 'N' ? 2 : 1) + chkstride1(A) + chkstride1(B) + chkstride1(C) ccall((@blasfunc($fname), libblas), Cvoid, (Ref{UInt8}, Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt}, Ref{$elty}, Ptr{$elty}, Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt}, Ref{$elty}, @@ -1369,12 +1391,12 @@ for (fname, elty) in ((:dsyr2k_,:Float64), end end end -function syr2k(uplo::Char, trans::Char, alpha::Number, A::StridedVecOrMat, B::StridedVecOrMat) +function syr2k(uplo::Char, trans::Char, alpha::Number, A::AbstractVecOrMat, B::AbstractVecOrMat) T = eltype(A) n = size(A, trans == 'N' ? 1 : 2) syr2k!(uplo, trans, convert(T,alpha), A, B, zero(T), similar(A, T, (n, n))) end -syr2k(uplo::Char, trans::Char, A::StridedVecOrMat, B::StridedVecOrMat) = syr2k(uplo, trans, one(eltype(A)), A, B) +syr2k(uplo::Char, trans::Char, A::AbstractVecOrMat, B::AbstractVecOrMat) = syr2k(uplo, trans, one(eltype(A)), A, B) for (fname, elty1, elty2) in ((:zher2k_,:ComplexF64,:Float64), (:cher2k_,:ComplexF32,:Float32)) @eval begin @@ -1389,11 +1411,14 @@ for (fname, elty1, elty2) in ((:zher2k_,:ComplexF64,:Float64), (:cher2k_,:Comple # .. Array Arguments .. # COMPLEX A(LDA,*),B(LDB,*),C(LDC,*) function her2k!(uplo::Char, trans::Char, alpha::($elty1), - A::StridedVecOrMat{$elty1}, B::StridedVecOrMat{$elty1}, - beta::($elty2), C::StridedMatrix{$elty1}) + A::AbstractVecOrMat{$elty1}, B::AbstractVecOrMat{$elty1}, + beta::($elty2), C::AbstractMatrix{$elty1}) n = checksquare(C) nn = size(A, trans == 'N' ? 1 : 2) if nn != n throw(DimensionMismatch("C has size ($n,$n), corresponding dimension of A is $nn")) end + chkstride1(A) + chkstride1(B) + chkstride1(C) k = size(A, trans == 'N' ? 2 : 1) ccall((@blasfunc($fname), libblas), Cvoid, (Ref{UInt8}, Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt}, @@ -1404,11 +1429,11 @@ for (fname, elty1, elty2) in ((:zher2k_,:ComplexF64,:Float64), (:cher2k_,:Comple beta, C, max(1,stride(C,2))) C end - function her2k(uplo::Char, trans::Char, alpha::($elty1), A::StridedVecOrMat{$elty1}, B::StridedVecOrMat{$elty1}) + function her2k(uplo::Char, trans::Char, alpha::($elty1), A::AbstractVecOrMat{$elty1}, B::AbstractVecOrMat{$elty1}) n = size(A, trans == 'N' ? 1 : 2) her2k!(uplo, trans, alpha, A, B, zero($elty2), similar(A, $elty1, (n,n))) end - her2k(uplo::Char, trans::Char, A::StridedVecOrMat{$elty1}, B::StridedVecOrMat{$elty1}) = her2k(uplo, trans, one($elty1), A, B) + her2k(uplo::Char, trans::Char, A::AbstractVecOrMat{$elty1}, B::AbstractVecOrMat{$elty1}) = her2k(uplo, trans, one($elty1), A, B) end end @@ -1474,12 +1499,14 @@ for (mmname, smname, elty) in # * .. Array Arguments .. # DOUBLE PRECISION A(LDA,*),B(LDB,*) function trmm!(side::Char, uplo::Char, transa::Char, diag::Char, alpha::Number, - A::StridedMatrix{$elty}, B::StridedMatrix{$elty}) + A::AbstractMatrix{$elty}, B::AbstractMatrix{$elty}) m, n = size(B) nA = checksquare(A) if nA != (side == 'L' ? m : n) throw(DimensionMismatch("size of A, $(size(A)), doesn't match $side size of B with dims, $(size(B))")) end + chkstride1(A) + chkstride1(B) ccall((@blasfunc($mmname), libblas), Cvoid, (Ref{UInt8}, Ref{UInt8}, Ref{UInt8}, Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt}, Ref{$elty}, Ptr{$elty}, Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt}), @@ -1488,7 +1515,7 @@ for (mmname, smname, elty) in B end function trmm(side::Char, uplo::Char, transa::Char, diag::Char, - alpha::$elty, A::StridedMatrix{$elty}, B::StridedMatrix{$elty}) + alpha::$elty, A::AbstractMatrix{$elty}, B::AbstractMatrix{$elty}) trmm!(side, uplo, transa, diag, alpha, A, copy(B)) end # SUBROUTINE DTRSM(SIDE,UPLO,TRANSA,DIAG,M,N,ALPHA,A,LDA,B,LDB) @@ -1499,12 +1526,14 @@ for (mmname, smname, elty) in # * .. Array Arguments .. # DOUBLE PRECISION A(LDA,*),B(LDB,*) function trsm!(side::Char, uplo::Char, transa::Char, diag::Char, - alpha::$elty, A::StridedMatrix{$elty}, B::StridedMatrix{$elty}) + alpha::$elty, A::AbstractMatrix{$elty}, B::AbstractMatrix{$elty}) m, n = size(B) k = checksquare(A) if k != (side == 'L' ? m : n) throw(DimensionMismatch("size of A is ($k,$k), size of B is ($m,$n), side is $side, and transa='$transa'")) end + chkstride1(A) + chkstride1(B) ccall((@blasfunc($smname), libblas), Cvoid, (Ref{UInt8}, Ref{UInt8}, Ref{UInt8}, Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt}, Ref{$elty}, Ptr{$elty}, @@ -1514,7 +1543,7 @@ for (mmname, smname, elty) in max(1,stride(A,2)), B, max(1,stride(B,2))) B end - function trsm(side::Char, uplo::Char, transa::Char, diag::Char, alpha::$elty, A::StridedMatrix{$elty}, B::StridedMatrix{$elty}) + function trsm(side::Char, uplo::Char, transa::Char, diag::Char, alpha::$elty, A::AbstractMatrix{$elty}, B::AbstractMatrix{$elty}) trsm!(side, uplo, transa, diag, alpha, A, copy(B)) end end diff --git a/base/linalg/dense.jl b/base/linalg/dense.jl index f27a901b02150..4a49b6b85ddc8 100644 --- a/base/linalg/dense.jl +++ b/base/linalg/dense.jl @@ -112,35 +112,21 @@ true isposdef(A::AbstractMatrix) = ishermitian(A) && isposdef(cholfact(Hermitian(A))) isposdef(x::Number) = imag(x)==0 && real(x) > 0 -""" - stride1(A) -> Int - -Return the distance between successive array elements -in dimension 1 in units of element size. - -# Examples -```jldoctest -julia> A = [1,2,3,4] -4-element Array{Int64,1}: - 1 - 2 - 3 - 4 - -julia> Base.LinAlg.stride1(A) -1 - -julia> B = view(A, 2:2:4) -2-element view(::Array{Int64,1}, 2:2:4) with eltype Int64: - 2 - 4 - -julia> Base.LinAlg.stride1(B) -2 -``` -""" -stride1(x::Array) = 1 -stride1(x::StridedVector) = stride(x, 1)::Int +# the definition of strides for Array{T,N} is tuple() if N = 0, otherwise it is +# a tuple containing 1 and a cumulative product of the first N-1 sizes +# this definition is also used for StridedReshapedArray and StridedReinterpretedArray +# which have the same memory storage as Array +function stride(a::Union{DenseArray,StridedReshapedArray,StridedReinterpretArray}, i::Int) + if i > ndims(a) + return length(a) + end + s = 1 + for n = 1:(i-1) + s *= size(a, n) + end + return s +end +strides(a::Union{DenseArray,StridedReshapedArray,StridedReinterpretArray}) = size_to_strides(1, size(a)...) function norm(x::StridedVector{T}, rx::Union{UnitRange{TI},AbstractRange{TI}}) where {T<:BlasFloat,TI<:Integer} if minimum(rx) < 1 || maximum(rx) > length(x) diff --git a/base/linalg/lapack.jl b/base/linalg/lapack.jl index dec2ed11b003f..27188a9b07858 100644 --- a/base/linalg/lapack.jl +++ b/base/linalg/lapack.jl @@ -84,7 +84,7 @@ end subsetrows(X::AbstractVector, Y::AbstractArray, k) = Y[1:k] subsetrows(X::AbstractMatrix, Y::AbstractArray, k) = Y[1:k, :] -function chkfinite(A::StridedMatrix) +function chkfinite(A::AbstractMatrix) for a in A if !isfinite(a) throw(ArgumentError("matrix contains Infs or NaNs")) @@ -117,7 +117,7 @@ for (gbtrf, gbtrs, elty) in # * .. Array Arguments .. # INTEGER IPIV( * ) # DOUBLE PRECISION AB( LDAB, * ) - function gbtrf!(kl::Integer, ku::Integer, m::Integer, AB::StridedMatrix{$elty}) + function gbtrf!(kl::Integer, ku::Integer, m::Integer, AB::AbstractMatrix{$elty}) chkstride1(AB) n = size(AB, 2) mnmn = min(m, n) @@ -139,8 +139,8 @@ for (gbtrf, gbtrs, elty) in # INTEGER IPIV( * ) # DOUBLE PRECISION AB( LDAB, * ), B( LDB, * ) function gbtrs!(trans::Char, kl::Integer, ku::Integer, m::Integer, - AB::StridedMatrix{$elty}, ipiv::StridedVector{BlasInt}, - B::StridedVecOrMat{$elty}) + AB::AbstractMatrix{$elty}, ipiv::AbstractVector{BlasInt}, + B::AbstractVecOrMat{$elty}) chkstride1(AB, B, ipiv) chktrans(trans) info = Ref{BlasInt}() @@ -168,7 +168,7 @@ subdiagonal containing a nonzero band, `ku` is the last superdiagonal containing one, and `m` is the first dimension of the matrix `AB`. Returns the LU factorization in-place and `ipiv`, the vector of pivots used. """ -gbtrf!(kl::Integer, ku::Integer, m::Integer, AB::StridedMatrix) +gbtrf!(kl::Integer, ku::Integer, m::Integer, AB::AbstractMatrix) """ gbtrs!(trans, kl, ku, m, AB, ipiv, B) @@ -179,7 +179,7 @@ first subdiagonal containing a nonzero band, `ku` is the last superdiagonal containing one, and `m` is the first dimension of the matrix `AB`. `ipiv` is the vector of pivots returned from `gbtrf!`. Returns the vector or matrix `X`, overwriting `B` in-place. """ -gbtrs!(trans::Char, kl::Integer, ku::Integer, m::Integer, AB::StridedMatrix, ipiv::StridedVector{BlasInt}, B::StridedVecOrMat) +gbtrs!(trans::Char, kl::Integer, ku::Integer, m::Integer, AB::AbstractMatrix, ipiv::AbstractVector{BlasInt}, B::AbstractVecOrMat) ## (GE) general matrices: balancing and back-transforming for (gebal, gebak, elty, relty) in @@ -194,7 +194,7 @@ for (gebal, gebak, elty, relty) in # INTEGER IHI, ILP, INFO, LDA, N # .. Array Arguments .. # DOUBLE PRECISION A( LDA, * ), SCALE( * ) - function gebal!(job::Char, A::StridedMatrix{$elty}) + function gebal!(job::Char, A::AbstractMatrix{$elty}) chkstride1(A) n = checksquare(A) chkfinite(A) # balancing routines don't support NaNs and Infs @@ -217,8 +217,8 @@ for (gebal, gebak, elty, relty) in # .. Array Arguments .. # DOUBLE PRECISION SCALE( * ), V( LDV, * ) function gebak!(job::Char, side::Char, - ilo::BlasInt, ihi::BlasInt, scale::StridedVector{$relty}, - V::StridedMatrix{$elty}) + ilo::BlasInt, ihi::BlasInt, scale::AbstractVector{$relty}, + V::AbstractMatrix{$elty}) chkstride1(scale, V) chkside(side) chkfinite(V) # balancing routines don't support NaNs and Infs @@ -244,7 +244,7 @@ and scaled). Modifies `A` in-place and returns `ilo`, `ihi`, and `scale`. If permuting was turned on, `A[i,j] = 0` if `j > i` and `1 < j < ilo` or `j > ihi`. `scale` contains information about the scaling/permutations performed. """ -gebal!(job::Char, A::StridedMatrix) +gebal!(job::Char, A::AbstractMatrix) """ gebak!(job, side, ilo, ihi, scale, V) @@ -254,7 +254,7 @@ the unscaled/unpermuted eigenvectors of the original matrix. Modifies `V` in-place. `side` can be `L` (left eigenvectors are transformed) or `R` (right eigenvectors are transformed). """ -gebak!(job::Char, side::Char, ilo::BlasInt, ihi::BlasInt, scale::StridedVector, V::StridedMatrix) +gebak!(job::Char, side::Char, ilo::BlasInt, ihi::BlasInt, scale::AbstractVector, V::AbstractMatrix) # (GE) general matrices, direct decompositions # @@ -276,7 +276,7 @@ for (gebrd, gelqf, geqlf, geqrf, geqp3, geqrt, geqrt3, gerqf, getrf, elty, relty # .. Array Arguments .. # DOUBLE PRECISION A( LDA, * ), D( * ), E( * ), TAUP( * ), # TAUQ( * ), WORK( * ) - function gebrd!(A::StridedMatrix{$elty}) + function gebrd!(A::AbstractMatrix{$elty}) chkstride1(A) m, n = size(A) k = min(m, n) @@ -309,7 +309,7 @@ for (gebrd, gelqf, geqlf, geqrf, geqp3, geqrt, geqrt3, gerqf, getrf, elty, relty # INTEGER INFO, LDA, LWORK, M, N # * .. Array Arguments .. # DOUBLE PRECISION A( LDA, * ), TAU( * ), WORK( * ) - function gelqf!(A::StridedMatrix{$elty}, tau::StridedVector{$elty}) + function gelqf!(A::AbstractMatrix{$elty}, tau::AbstractVector{$elty}) chkstride1(A,tau) m = BlasInt(size(A, 1)) n = BlasInt(size(A, 2)) @@ -339,7 +339,7 @@ for (gebrd, gelqf, geqlf, geqrf, geqp3, geqrt, geqrt3, gerqf, getrf, elty, relty # INTEGER INFO, LDA, LWORK, M, N # * .. Array Arguments .. # DOUBLE PRECISION A( LDA, * ), TAU( * ), WORK( * ) - function geqlf!(A::StridedMatrix{$elty}, tau::StridedVector{$elty}) + function geqlf!(A::AbstractMatrix{$elty}, tau::AbstractVector{$elty}) chkstride1(A,tau) m = BlasInt(size(A, 1)) n = BlasInt(size(A, 2)) @@ -370,7 +370,7 @@ for (gebrd, gelqf, geqlf, geqrf, geqp3, geqrt, geqrt3, gerqf, getrf, elty, relty # * .. Array Arguments .. # INTEGER JPVT( * ) # DOUBLE PRECISION A( LDA, * ), TAU( * ), WORK( * ) - function geqp3!(A::StridedMatrix{$elty}, jpvt::StridedVector{BlasInt}, tau::StridedVector{$elty}) + function geqp3!(A::AbstractMatrix{$elty}, jpvt::AbstractVector{BlasInt}, tau::AbstractVector{$elty}) chkstride1(A,jpvt,tau) m,n = size(A) if length(tau) != min(m,n) @@ -417,7 +417,7 @@ for (gebrd, gelqf, geqlf, geqrf, geqp3, geqrt, geqrt3, gerqf, getrf, elty, relty return A, tau, jpvt end - function geqrt!(A::StridedMatrix{$elty}, T::StridedMatrix{$elty}) + function geqrt!(A::AbstractMatrix{$elty}, T::AbstractMatrix{$elty}) chkstride1(A) m, n = size(A) minmn = min(m, n) @@ -441,7 +441,7 @@ for (gebrd, gelqf, geqlf, geqrf, geqp3, geqrt, geqrt3, gerqf, getrf, elty, relty A, T end - function geqrt3!(A::StridedMatrix{$elty}, T::StridedMatrix{$elty}) + function geqrt3!(A::AbstractMatrix{$elty}, T::AbstractMatrix{$elty}) chkstride1(A) chkstride1(T) m, n = size(A) @@ -470,7 +470,7 @@ for (gebrd, gelqf, geqlf, geqrf, geqp3, geqrt, geqrt3, gerqf, getrf, elty, relty # INTEGER INFO, LDA, LWORK, M, N # * .. Array Arguments .. # DOUBLE PRECISION A( LDA, * ), TAU( * ), WORK( * ) - function geqrf!(A::StridedMatrix{$elty}, tau::StridedVector{$elty}) + function geqrf!(A::AbstractMatrix{$elty}, tau::AbstractVector{$elty}) chkstride1(A,tau) m, n = size(A) if length(tau) != min(m,n) @@ -498,7 +498,7 @@ for (gebrd, gelqf, geqlf, geqrf, geqp3, geqrt, geqrt3, gerqf, getrf, elty, relty # INTEGER INFO, LDA, LWORK, M, N # * .. Array Arguments .. # DOUBLE PRECISION A( LDA, * ), TAU( * ), WORK( * ) - function gerqf!(A::StridedMatrix{$elty},tau::StridedVector{$elty}) + function gerqf!(A::AbstractMatrix{$elty},tau::AbstractVector{$elty}) chkstride1(A,tau) m, n = size(A) if length(tau) != min(m,n) @@ -527,7 +527,7 @@ for (gebrd, gelqf, geqlf, geqrf, geqp3, geqrt, geqrt3, gerqf, getrf, elty, relty # * .. Array Arguments .. # INTEGER IPIV( * ) # DOUBLE PRECISION A( LDA, * ) - function getrf!(A::StridedMatrix{$elty}) + function getrf!(A::AbstractMatrix{$elty}) chkstride1(A) m, n = size(A) lda = max(1,stride(A, 2)) @@ -552,7 +552,7 @@ containing the off-diagonal elements of `B`; `tauq`, containing the elementary reflectors representing `Q`; and `taup`, containing the elementary reflectors representing `P`. """ -gebrd!(A::StridedMatrix) +gebrd!(A::AbstractMatrix) """ gelqf!(A, tau) @@ -564,7 +564,7 @@ must have length greater than or equal to the smallest dimension of `A`. Returns `A` and `tau` modified in-place. """ -gelqf!(A::StridedMatrix, tau::StridedVector) +gelqf!(A::AbstractMatrix, tau::AbstractVector) """ geqlf!(A, tau) @@ -575,7 +575,7 @@ must have length greater than or equal to the smallest dimension of `A`. Returns `A` and `tau` modified in-place. """ -geqlf!(A::StridedMatrix, tau::StridedVector) +geqlf!(A::AbstractMatrix, tau::AbstractVector) """ geqp3!(A, jpvt, tau) @@ -588,7 +588,7 @@ smallest dimension of `A`. `A`, `jpvt`, and `tau` are modified in-place. """ -geqp3!(A::StridedMatrix, jpvt::StridedVector{BlasInt}, tau::StridedVector) +geqp3!(A::AbstractMatrix, jpvt::AbstractVector{BlasInt}, tau::AbstractVector) """ geqrt!(A, T) @@ -601,7 +601,7 @@ dimension of `A`. Returns `A` and `T` modified in-place. """ -geqrt!(A::StridedMatrix, T::StridedMatrix) +geqrt!(A::AbstractMatrix, T::AbstractMatrix) """ geqrt3!(A, T) @@ -614,7 +614,7 @@ equal the smallest dimension of `A`. Returns `A` and `T` modified in-place. """ -geqrt3!(A::StridedMatrix, T::StridedMatrix) +geqrt3!(A::AbstractMatrix, T::AbstractMatrix) """ geqrf!(A, tau) @@ -625,7 +625,7 @@ must have length greater than or equal to the smallest dimension of `A`. Returns `A` and `tau` modified in-place. """ -geqrf!(A::StridedMatrix, tau::StridedVector) +geqrf!(A::AbstractMatrix, tau::AbstractVector) """ gerqf!(A, tau) @@ -636,7 +636,7 @@ must have length greater than or equal to the smallest dimension of `A`. Returns `A` and `tau` modified in-place. """ -gerqf!(A::StridedMatrix, tau::StridedVector) +gerqf!(A::AbstractMatrix, tau::AbstractVector) """ getrf!(A) -> (A, ipiv, info) @@ -647,7 +647,7 @@ Returns `A`, modified in-place, `ipiv`, the pivoting information, and an `info` code which indicates success (`info = 0`), a singular value in `U` (`info = i`, in which case `U[i,i]` is singular), or an error code (`info < 0`). """ -getrf!(A::StridedMatrix, tau::StridedVector) +getrf!(A::AbstractMatrix, tau::AbstractVector) """ gelqf!(A) -> (A, tau) @@ -657,7 +657,7 @@ Compute the `LQ` factorization of `A`, `A = LQ`. Returns `A`, modified in-place, and `tau`, which contains scalars which parameterize the elementary reflectors of the factorization. """ -gelqf!(A::StridedMatrix{<:BlasFloat}) = ((m,n) = size(A); gelqf!(A, similar(A, min(m, n)))) +gelqf!(A::AbstractMatrix{<:BlasFloat}) = ((m,n) = size(A); gelqf!(A, similar(A, min(m, n)))) """ geqlf!(A) -> (A, tau) @@ -667,7 +667,7 @@ Compute the `QL` factorization of `A`, `A = QL`. Returns `A`, modified in-place, and `tau`, which contains scalars which parameterize the elementary reflectors of the factorization. """ -geqlf!(A::StridedMatrix{<:BlasFloat}) = ((m,n) = size(A); geqlf!(A, similar(A, min(m, n)))) +geqlf!(A::AbstractMatrix{<:BlasFloat}) = ((m,n) = size(A); geqlf!(A, similar(A, min(m, n)))) """ geqrt!(A, nb) -> (A, T) @@ -679,7 +679,7 @@ Returns `A`, modified in-place, and `T`, which contains upper triangular block reflectors which parameterize the elementary reflectors of the factorization. """ -geqrt!(A::StridedMatrix{<:BlasFloat}, nb::Integer) = geqrt!(A, similar(A, nb, minimum(size(A)))) +geqrt!(A::AbstractMatrix{<:BlasFloat}, nb::Integer) = geqrt!(A, similar(A, nb, minimum(size(A)))) """ geqrt3!(A) -> (A, T) @@ -689,7 +689,7 @@ Recursively computes the blocked `QR` factorization of `A`, `A = QR`. Returns `A`, modified in-place, and `T`, which contains upper triangular block reflectors which parameterize the elementary reflectors of the factorization. """ -geqrt3!(A::StridedMatrix{<:BlasFloat}) = (n = size(A, 2); geqrt3!(A, similar(A, n, n))) +geqrt3!(A::AbstractMatrix{<:BlasFloat}) = (n = size(A, 2); geqrt3!(A, similar(A, n, n))) """ geqrf!(A) -> (A, tau) @@ -699,7 +699,7 @@ Compute the `QR` factorization of `A`, `A = QR`. Returns `A`, modified in-place, and `tau`, which contains scalars which parameterize the elementary reflectors of the factorization. """ -geqrf!(A::StridedMatrix{<:BlasFloat}) = ((m,n) = size(A); geqrf!(A, similar(A, min(m, n)))) +geqrf!(A::AbstractMatrix{<:BlasFloat}) = ((m,n) = size(A); geqrf!(A, similar(A, min(m, n)))) """ gerqf!(A) -> (A, tau) @@ -709,7 +709,7 @@ Compute the `RQ` factorization of `A`, `A = RQ`. Returns `A`, modified in-place, and `tau`, which contains scalars which parameterize the elementary reflectors of the factorization. """ -gerqf!(A::StridedMatrix{<:BlasFloat}) = ((m,n) = size(A); gerqf!(A, similar(A, min(m, n)))) +gerqf!(A::AbstractMatrix{<:BlasFloat}) = ((m,n) = size(A); gerqf!(A, similar(A, min(m, n)))) """ geqp3!(A, jpvt) -> (A, jpvt, tau) @@ -721,7 +721,7 @@ greater than or equal to `n` if `A` is an `(m x n)` matrix. Returns `A` and `jpvt`, modified in-place, and `tau`, which stores the elementary reflectors. """ -function geqp3!(A::StridedMatrix{<:BlasFloat}, jpvt::StridedVector{BlasInt}) +function geqp3!(A::AbstractMatrix{<:BlasFloat}, jpvt::AbstractVector{BlasInt}) m, n = size(A) geqp3!(A, jpvt, similar(A, min(m, n))) end @@ -734,7 +734,7 @@ Compute the pivoted `QR` factorization of `A`, `AP = QR` using BLAS level 3. Returns `A`, modified in-place, `jpvt`, which represents the pivoting matrix `P`, and `tau`, which stores the elementary reflectors. """ -function geqp3!(A::StridedMatrix{<:BlasFloat}) +function geqp3!(A::AbstractMatrix{<:BlasFloat}) m, n = size(A) geqp3!(A, zeros(BlasInt, n), similar(A, min(m, n))) end @@ -753,12 +753,13 @@ for (tzrzf, ormrz, elty) in # .. # .. Array Arguments .. # COMPLEX*16 A( LDA, * ), TAU( * ), WORK( * ) - function tzrzf!(A::StridedMatrix{$elty}) + function tzrzf!(A::AbstractMatrix{$elty}) + chkstride1(A) m, n = size(A) if n < m throw(DimensionMismatch("input matrix A has dimensions ($m,$n), but cannot have fewer columns than rows")) end - lda = max(1, m) + lda = max(1, stride(A,2)) tau = similar(A, $elty, m) work = Vector{$elty}(uninitialized, 1) lwork = BlasInt(-1) @@ -787,11 +788,11 @@ for (tzrzf, ormrz, elty) in # .. # .. Array Arguments .. # COMPLEX*16 A( LDA, * ), C( LDC, * ), TAU( * ), WORK( * ) - function ormrz!(side::Char, trans::Char, A::StridedMatrix{$elty}, - tau::StridedVector{$elty}, C::StridedMatrix{$elty}) + function ormrz!(side::Char, trans::Char, A::AbstractMatrix{$elty}, + tau::AbstractVector{$elty}, C::AbstractMatrix{$elty}) chktrans(trans) chkside(side) - chkstride1(tau) + chkstride1(A, tau, C) m, n = size(C) k = length(tau) l = size(A, 2) - size(A, 1) @@ -831,7 +832,7 @@ can be unmodified (`trans = N`), transposed (`trans = T`), or conjugate transposed (`trans = C`). Returns matrix `C` which is modified in-place with the result of the multiplication. """ -ormrz!(side::Char, trans::Char, A::StridedMatrix, tau::StridedVector, C::StridedMatrix) +ormrz!(side::Char, trans::Char, A::AbstractMatrix, tau::AbstractVector, C::AbstractMatrix) """ tzrzf!(A) -> (A, tau) @@ -840,7 +841,7 @@ Transforms the upper trapezoidal matrix `A` to upper triangular form in-place. Returns `A` and `tau`, the scalar parameters for the elementary reflectors of the transformation. """ -tzrzf!(A::StridedMatrix) +tzrzf!(A::AbstractMatrix) ## (GE) general matrices, solvers with factorization, solver and inverse for (gels, gesv, getrs, getri, elty) in @@ -853,7 +854,7 @@ for (gels, gesv, getrs, getri, elty) in # * .. Scalar Arguments .. # CHARACTER TRANS # INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS - function gels!(trans::Char, A::StridedMatrix{$elty}, B::StridedVecOrMat{$elty}) + function gels!(trans::Char, A::AbstractMatrix{$elty}, B::AbstractVecOrMat{$elty}) chktrans(trans) chkstride1(A, B) btrn = trans == 'T' @@ -897,7 +898,7 @@ for (gels, gesv, getrs, getri, elty) in # * .. Array Arguments .. # INTEGER IPIV( * ) # DOUBLE PRECISION A( LDA, * ), B( LDB, * ) - function gesv!(A::StridedMatrix{$elty}, B::StridedVecOrMat{$elty}) + function gesv!(A::AbstractMatrix{$elty}, B::AbstractVecOrMat{$elty}) chkstride1(A, B) n = checksquare(A) if size(B,1) != n @@ -920,7 +921,7 @@ for (gels, gesv, getrs, getri, elty) in # .. Array Arguments .. # INTEGER IPIV( * ) # DOUBLE PRECISION A( LDA, * ), B( LDB, * ) - function getrs!(trans::Char, A::StridedMatrix{$elty}, ipiv::StridedVector{BlasInt}, B::StridedVecOrMat{$elty}) + function getrs!(trans::Char, A::AbstractMatrix{$elty}, ipiv::AbstractVector{BlasInt}, B::AbstractVecOrMat{$elty}) chktrans(trans) chkstride1(A, B, ipiv) n = checksquare(A) @@ -943,7 +944,7 @@ for (gels, gesv, getrs, getri, elty) in #* .. Array Arguments .. # INTEGER IPIV( * ) # DOUBLE PRECISION A( LDA, * ), WORK( * ) - function getri!(A::StridedMatrix{$elty}, ipiv::StridedVector{BlasInt}) + function getri!(A::AbstractMatrix{$elty}, ipiv::AbstractVector{BlasInt}) chkstride1(A, ipiv) n = checksquare(A) if n != length(ipiv) @@ -979,7 +980,7 @@ may be one of `N` (no modification), `T` (transpose), or `C` (conjugate transpose). `gels!` searches for the minimum norm/least squares solution. `A` may be under or over determined. The solution is returned in `B`. """ -gels!(trans::Char, A::StridedMatrix, B::StridedVecOrMat) +gels!(trans::Char, A::AbstractMatrix, B::AbstractVecOrMat) """ gesv!(A, B) -> (B, A, ipiv) @@ -989,7 +990,7 @@ the `LU` factorization of `A`. `A` is overwritten with its `LU` factorization and `B` is overwritten with the solution `X`. `ipiv` contains the pivoting information for the `LU` factorization of `A`. """ -gesv!(A::StridedMatrix, B::StridedVecOrMat) +gesv!(A::AbstractMatrix, B::AbstractVecOrMat) """ getrs!(trans, A, ipiv, B) @@ -1000,7 +1001,7 @@ is the `LU` factorization from `getrf!`, with `ipiv` the pivoting information. `trans` may be one of `N` (no modification), `T` (transpose), or `C` (conjugate transpose). """ -getrs!(trans::Char, A::StridedMatrix, ipiv::StridedVector{BlasInt}, B::StridedVecOrMat) +getrs!(trans::Char, A::AbstractMatrix, ipiv::AbstractVector{BlasInt}, B::AbstractVecOrMat) """ getri!(A, ipiv) @@ -1010,7 +1011,7 @@ Computes the inverse of `A`, using its `LU` factorization found by contains the `LU` factorization of `getrf!`. `A` is overwritten with its inverse. """ -getri!(A::StridedMatrix, ipiv::StridedVector{BlasInt}) +getri!(A::AbstractMatrix, ipiv::AbstractVector{BlasInt}) for (gesvx, elty) in ((:dgesvx_,:Float64), @@ -1031,11 +1032,11 @@ for (gesvx, elty) in # $ BERR( * ), C( * ), FERR( * ), R( * ), # $ WORK( * ), X( LDX, * # - function gesvx!(fact::Char, trans::Char, A::StridedMatrix{$elty}, - AF::StridedMatrix{$elty}, ipiv::StridedVector{BlasInt}, equed::Char, - R::StridedVector{$elty}, C::StridedVector{$elty}, B::StridedVecOrMat{$elty}) + function gesvx!(fact::Char, trans::Char, A::AbstractMatrix{$elty}, + AF::AbstractMatrix{$elty}, ipiv::AbstractVector{BlasInt}, equed::Char, + R::AbstractVector{$elty}, C::AbstractVector{$elty}, B::AbstractVecOrMat{$elty}) chktrans(trans) - chkstride1(ipiv, R, C) + chkstride1(ipiv, R, C, B) n = checksquare(A) lda = stride(A,2) n = checksquare(AF) @@ -1067,7 +1068,7 @@ for (gesvx, elty) in X, equed, R, C, B, rcond[], ferr, berr, work[1] end - function gesvx!(A::StridedMatrix{$elty}, B::StridedVecOrMat{$elty}) + function gesvx!(A::AbstractMatrix{$elty}, B::AbstractVecOrMat{$elty}) n = size(A,1) X, equed, R, C, B, rcond, ferr, berr, rpgf = gesvx!('N', 'N', A, @@ -1100,11 +1101,11 @@ for (gesvx, elty, relty) in # $ RWORK( * ) # COMPLEX*16 A( LDA, * ), AF( LDAF, * ), B( LDB, * ), # $ WORK( * ), X( LDX, * ) - function gesvx!(fact::Char, trans::Char, A::StridedMatrix{$elty}, - AF::StridedMatrix{$elty}, ipiv::StridedVector{BlasInt}, equed::Char, - R::StridedVector{$relty}, C::StridedVector{$relty}, B::StridedVecOrMat{$elty}) + function gesvx!(fact::Char, trans::Char, A::AbstractMatrix{$elty}, + AF::AbstractMatrix{$elty}, ipiv::AbstractVector{BlasInt}, equed::Char, + R::AbstractVector{$relty}, C::AbstractVector{$relty}, B::AbstractVecOrMat{$elty}) chktrans(trans) - chkstride1(ipiv, R, C) + chkstride1(A, AF, ipiv, R, C, B) n = checksquare(A) lda = stride(A,2) n = checksquare(AF) @@ -1137,7 +1138,7 @@ for (gesvx, elty, relty) in end #Wrapper for the no-equilibration, no-transpose calculation - function gesvx!(A::StridedMatrix{$elty}, B::StridedVecOrMat{$elty}) + function gesvx!(A::AbstractMatrix{$elty}, B::AbstractVecOrMat{$elty}) n = size(A,1) X, equed, R, C, B, rcond, ferr, berr, rpgf = gesvx!('N', 'N', A, @@ -1176,15 +1177,15 @@ condition number of `A` after equilbrating; `ferr`, the forward error bound for each solution vector in `X`; `berr`, the forward error bound for each solution vector in `X`; and `work`, the reciprocal pivot growth factor. """ -gesvx!(fact::Char, trans::Char, A::StridedMatrix, AF::StridedMatrix, - ipiv::StridedVector{BlasInt}, equed::Char, R::StridedVector, C::StridedVector, B::StridedVecOrMat) +gesvx!(fact::Char, trans::Char, A::AbstractMatrix, AF::AbstractMatrix, + ipiv::AbstractVector{BlasInt}, equed::Char, R::AbstractVector, C::AbstractVector, B::AbstractVecOrMat) """ gesvx!(A, B) The no-equilibration, no-transpose simplification of `gesvx!`. """ -gesvx!(A::StridedMatrix, B::StridedVecOrMat) +gesvx!(A::AbstractMatrix, B::AbstractVecOrMat) for (gelsd, gelsy, elty) in ((:dgelsd_,:dgelsy_,:Float64), @@ -1199,7 +1200,7 @@ for (gelsd, gelsy, elty) in # * .. Array Arguments .. # INTEGER IWORK( * ) # DOUBLE PRECISION A( LDA, * ), B( LDB, * ), S( * ), WORK( * ) - function gelsd!(A::StridedMatrix{$elty}, B::StridedVecOrMat{$elty}, rcond::Real=-one($elty)) + function gelsd!(A::AbstractMatrix{$elty}, B::AbstractVecOrMat{$elty}, rcond::Real=-one($elty)) chkstride1(A, B) m, n = size(A) if size(B, 1) != m @@ -1241,8 +1242,8 @@ for (gelsd, gelsy, elty) in # * .. Array Arguments .. # INTEGER JPVT( * ) # DOUBLE PRECISION A( LDA, * ), B( LDB, * ), WORK( * ) - function gelsy!(A::StridedMatrix{$elty}, B::StridedVecOrMat{$elty}, rcond::Real=eps($elty)) - chkstride1(A, B) + function gelsy!(A::AbstractMatrix{$elty}, B::AbstractVecOrMat{$elty}, rcond::Real=eps($elty)) + chkstride1(A) m = size(A, 1) n = size(A, 2) nrhs = size(B, 2) @@ -1250,8 +1251,8 @@ for (gelsd, gelsy, elty) in throw(DimensionMismatch("B has leading dimension $(size(B,1)) but needs $m")) end newB = [B; zeros($elty, max(0, n - size(B, 1)), size(B, 2))] - lda = max(1, m) - ldb = max(1, m, n) + lda = max(1, stride(A,2)) + ldb = max(1, stride(newB,2)) jpvt = zeros(BlasInt, n) rnk = Ref{BlasInt}() work = Vector{$elty}(uninitialized, 1) @@ -1292,7 +1293,7 @@ for (gelsd, gelsy, elty, relty) in # INTEGER IWORK( * ) # DOUBLE PRECISION RWORK( * ), S( * ) # COMPLEX*16 A( LDA, * ), B( LDB, * ), WORK( * ) - function gelsd!(A::StridedMatrix{$elty}, B::StridedVecOrMat{$elty}, rcond::Real=-one($relty)) + function gelsd!(A::AbstractMatrix{$elty}, B::AbstractVecOrMat{$elty}, rcond::Real=-one($relty)) chkstride1(A, B) m, n = size(A) if size(B, 1) != m @@ -1337,7 +1338,7 @@ for (gelsd, gelsy, elty, relty) in # INTEGER JPVT( * ) # DOUBLE PRECISION RWORK( * ) # COMPLEX*16 A( LDA, * ), B( LDB, * ), WORK( * ) - function gelsy!(A::StridedMatrix{$elty}, B::StridedVecOrMat{$elty}, rcond::Real=eps($relty)) + function gelsy!(A::AbstractMatrix{$elty}, B::AbstractVecOrMat{$elty}, rcond::Real=eps($relty)) chkstride1(A, B) m, n = size(A) nrhs = size(B, 2) @@ -1383,7 +1384,7 @@ is overwritten with the solution `X`. Singular values below `rcond` will be treated as zero. Returns the solution in `B` and the effective rank of `A` in `rnk`. """ -gelsd!(A::StridedMatrix, B::StridedVecOrMat, rcond::Real) +gelsd!(A::AbstractMatrix, B::AbstractVecOrMat, rcond::Real) """ gelsy!(A, B, rcond) -> (B, rnk) @@ -1394,7 +1395,7 @@ is overwritten with the solution `X`. Singular values below `rcond` will be treated as zero. Returns the solution in `B` and the effective rank of `A` in `rnk`. """ -gelsy!(A::StridedMatrix, B::StridedVecOrMat, rcond::Real) +gelsy!(A::AbstractMatrix, B::AbstractVecOrMat, rcond::Real) for (gglse, elty) in ((:dgglse_, :Float64), (:sgglse_, :Float32), @@ -1409,9 +1410,9 @@ for (gglse, elty) in ((:dgglse_, :Float64), # * .. Array Arguments .. # DOUBLE PRECISION A( LDA, * ), B( LDB, * ), C( * ), D( * ), # $ WORK( * ), X( * ) - function gglse!(A::StridedMatrix{$elty}, c::StridedVector{$elty}, - B::StridedMatrix{$elty}, d::StridedVector{$elty}) - chkstride1(A, B) + function gglse!(A::AbstractMatrix{$elty}, c::AbstractVector{$elty}, + B::AbstractMatrix{$elty}, d::AbstractVector{$elty}) + chkstride1(A, c, B, d) m, n = size(A) p = size(B, 1) if size(B, 2) != n @@ -1453,7 +1454,7 @@ Solves the equation `A * x = c` where `x` is subject to the equality constraint `B * x = d`. Uses the formula `||c - A*x||^2 = 0` to solve. Returns `X` and the residual sum-of-squares. """ -gglse!(A::StridedMatrix, c::StridedVector, B::StridedMatrix, d::StridedVector) +gglse!(A::AbstractMatrix, c::AbstractVector, B::AbstractMatrix, d::AbstractVector) # (GE) general matrices eigenvalue-eigenvector and singular value decompositions for (geev, gesvd, gesdd, ggsvd, elty, relty) in @@ -1470,7 +1471,7 @@ for (geev, gesvd, gesdd, ggsvd, elty, relty) in # * .. Array Arguments .. # DOUBLE PRECISION A( LDA, * ), VL( LDVL, * ), VR( LDVR, * ), # $ WI( * ), WORK( * ), WR( * ) - function geev!(jobvl::Char, jobvr::Char, A::StridedMatrix{$elty}) + function geev!(jobvl::Char, jobvr::Char, A::AbstractMatrix{$elty}) chkstride1(A) n = checksquare(A) chkfinite(A) # balancing routines don't support NaNs and Infs @@ -1526,7 +1527,7 @@ for (geev, gesvd, gesdd, ggsvd, elty, relty) in # INTEGER IWORK( * ) # DOUBLE PRECISION A( LDA, * ), S( * ), U( LDU, * ), # VT( LDVT, * ), WORK( * ) - function gesdd!(job::Char, A::StridedMatrix{$elty}) + function gesdd!(job::Char, A::AbstractMatrix{$elty}) chkstride1(A) m, n = size(A) minmn = min(m, n) @@ -1606,7 +1607,7 @@ for (geev, gesvd, gesdd, ggsvd, elty, relty) in # * .. Array Arguments .. # DOUBLE PRECISION A( LDA, * ), S( * ), U( LDU, * ), # $ VT( LDVT, * ), WORK( * ) - function gesvd!(jobu::Char, jobvt::Char, A::StridedMatrix{$elty}) + function gesvd!(jobu::Char, jobvt::Char, A::AbstractMatrix{$elty}) chkstride1(A) m, n = size(A) minmn = min(m, n) @@ -1674,7 +1675,7 @@ for (geev, gesvd, gesdd, ggsvd, elty, relty) in # DOUBLE PRECISION ALPHA( * ), BETA( * ), RWORK( * ) # COMPLEX*16 A( LDA, * ), B( LDB, * ), Q( LDQ, * ), # $ U( LDU, * ), V( LDV, * ), WORK( * ) - function ggsvd!(jobu::Char, jobv::Char, jobq::Char, A::StridedMatrix{$elty}, B::StridedMatrix{$elty}) + function ggsvd!(jobu::Char, jobv::Char, jobq::Char, A::AbstractMatrix{$elty}, B::AbstractMatrix{$elty}) chkstride1(A, B) m, n = size(A) if size(B, 2) != n @@ -1749,7 +1750,7 @@ aren't computed. If `jobvl = V` or `jobvr = V`, the corresponding eigenvectors are computed. Returns the eigenvalues in `W`, the right eigenvectors in `VR`, and the left eigenvectors in `VL`. """ -geev!(jobvl::Char, jobvr::Char, A::StridedMatrix) +geev!(jobvl::Char, jobvr::Char, A::AbstractMatrix) """ gesdd!(job, A) -> (U, S, VT) @@ -1761,7 +1762,7 @@ are computed. If `job = O`, `A` is overwritten with the columns of (thin) `U` and the rows of (thin) `V'`. If `job = S`, the columns of (thin) `U` and the rows of (thin) `V'` are computed and returned separately. """ -gesdd!(job::Char, A::StridedMatrix) +gesdd!(job::Char, A::AbstractMatrix) """ gesvd!(jobu, jobvt, A) -> (U, S, VT) @@ -1777,7 +1778,7 @@ computed and returned separately. `jobu` and `jobvt` can't both be `O`. Returns `U`, `S`, and `Vt`, where `S` are the singular values of `A`. """ -gesvd!(jobu::Char, jobvt::Char, A::StridedMatrix) +gesvd!(jobu::Char, jobvt::Char, A::AbstractMatrix) """ ggsvd!(jobu, jobv, jobq, A, B) -> (U, V, Q, alpha, beta, k, l, R) @@ -1790,13 +1791,13 @@ the orthogonal/unitary matrix `Q` is computed. If `jobu`, `jobv` or `jobq` is `N`, that matrix is not computed. This function is only available in LAPACK versions prior to 3.6.0. """ -ggsvd!(jobu::Char, jobv::Char, jobq::Char, A::StridedMatrix, B::StridedMatrix) +ggsvd!(jobu::Char, jobv::Char, jobq::Char, A::AbstractMatrix, B::AbstractMatrix) for (f, elty) in ((:dggsvd3_, :Float64), (:sggsvd3_, :Float32)) @eval begin - function ggsvd3!(jobu::Char, jobv::Char, jobq::Char, A::StridedMatrix{$elty}, B::StridedMatrix{$elty}) + function ggsvd3!(jobu::Char, jobv::Char, jobq::Char, A::AbstractMatrix{$elty}, B::AbstractMatrix{$elty}) chkstride1(A, B) m, n = size(A) if size(B, 2) != n @@ -1852,7 +1853,7 @@ end for (f, elty, relty) in ((:zggsvd3_, :ComplexF64, :Float64), (:cggsvd3_, :ComplexF32, :Float32)) @eval begin - function ggsvd3!(jobu::Char, jobv::Char, jobq::Char, A::StridedMatrix{$elty}, B::StridedMatrix{$elty}) + function ggsvd3!(jobu::Char, jobv::Char, jobq::Char, A::AbstractMatrix{$elty}, B::AbstractMatrix{$elty}) chkstride1(A, B) m, n = size(A) if size(B, 2) != n @@ -1939,7 +1940,7 @@ for (geevx, ggev, elty) in # DOUBLE PRECISION A( LDA, * ), RCONDE( * ), RCONDV( * ), # $ SCALE( * ), VL( LDVL, * ), VR( LDVR, * ), # $ WI( * ), WORK( * ), WR( * ) - function geevx!(balanc::Char, jobvl::Char, jobvr::Char, sense::Char, A::StridedMatrix{$elty}) + function geevx!(balanc::Char, jobvl::Char, jobvr::Char, sense::Char, A::AbstractMatrix{$elty}) n = checksquare(A) chkfinite(A) # balancing routines don't support NaNs and Infs lda = max(1,stride(A,2)) @@ -2017,7 +2018,7 @@ for (geevx, ggev, elty) in # DOUBLE PRECISION A( LDA, * ), ALPHAI( * ), ALPHAR( * ), # $ B( LDB, * ), BETA( * ), VL( LDVL, * ), # $ VR( LDVR, * ), WORK( * ) - function ggev!(jobvl::Char, jobvr::Char, A::StridedMatrix{$elty}, B::StridedMatrix{$elty}) + function ggev!(jobvl::Char, jobvr::Char, A::AbstractMatrix{$elty}, B::AbstractMatrix{$elty}) chkstride1(A,B) n, m = checksquare(A,B) if n != m @@ -2090,7 +2091,7 @@ for (geevx, ggev, elty, relty) in # $ SCALE( * ) # COMPLEX*16 A( LDA, * ), VL( LDVL, * ), VR( LDVR, * ), # $ W( * ), WORK( * ) - function geevx!(balanc::Char, jobvl::Char, jobvr::Char, sense::Char, A::StridedMatrix{$elty}) + function geevx!(balanc::Char, jobvl::Char, jobvr::Char, sense::Char, A::AbstractMatrix{$elty}) n = checksquare(A) chkfinite(A) # balancing routines don't support NaNs and Infs lda = max(1,stride(A,2)) @@ -2163,7 +2164,7 @@ for (geevx, ggev, elty, relty) in # COMPLEX*16 A( LDA, * ), ALPHA( * ), B( LDB, * ), # $ BETA( * ), VL( LDVL, * ), VR( LDVR, * ), # $ WORK( * ) - function ggev!(jobvl::Char, jobvr::Char, A::StridedMatrix{$elty}, B::StridedMatrix{$elty}) + function ggev!(jobvl::Char, jobvr::Char, A::AbstractMatrix{$elty}, B::AbstractMatrix{$elty}) chkstride1(A, B) n, m = checksquare(A, B) if n != m @@ -2235,7 +2236,7 @@ condition numbers are computed for the right eigenvectors and the eigenvectors. If `sense = E,B`, the right and left eigenvectors must be computed. """ -geevx!(balanc::Char, jobvl::Char, jobvr::Char, sense::Char, A::StridedMatrix) +geevx!(balanc::Char, jobvl::Char, jobvr::Char, sense::Char, A::AbstractMatrix) """ ggev!(jobvl, jobvr, A, B) -> (alpha, beta, vl, vr) @@ -2245,7 +2246,7 @@ the left eigenvectors aren't computed. If `jobvr = N`, the right eigenvectors aren't computed. If `jobvl = V` or `jobvr = V`, the corresponding eigenvectors are computed. """ -ggev!(jobvl::Char, jobvr::Char, A::StridedMatrix, B::StridedMatrix) +ggev!(jobvl::Char, jobvr::Char, A::AbstractMatrix, B::AbstractMatrix) # One step incremental condition estimation of max/min singular values for (laic1, elty) in @@ -2260,8 +2261,8 @@ for (laic1, elty) in # .. # .. Array Arguments .. # DOUBLE PRECISION W( J ), X( J ) - function laic1!(job::Integer, x::StridedVector{$elty}, - sest::$elty, w::StridedVector{$elty}, gamma::$elty) + function laic1!(job::Integer, x::AbstractVector{$elty}, + sest::$elty, w::AbstractVector{$elty}, gamma::$elty) j = length(x) if j != length(w) throw(DimensionMismatch("vectors must have same length, but length of x is $j and length of w is $(length(w))")) @@ -2293,8 +2294,8 @@ for (laic1, elty, relty) in # .. # .. Array Arguments .. # COMPLEX*16 W( J ), X( J ) - function laic1!(job::Integer, x::StridedVector{$elty}, - sest::$relty, w::StridedVector{$elty}, gamma::$elty) + function laic1!(job::Integer, x::AbstractVector{$elty}, + sest::$relty, w::AbstractVector{$elty}, gamma::$elty) j = length(x) if j != length(w) throw(DimensionMismatch("vectors must have same length, but length of x is $j and length of w is $(length(w))")) @@ -2326,8 +2327,8 @@ for (gtsv, gttrf, gttrs, elty) in # INTEGER INFO, LDB, N, NRHS # .. Array Arguments .. # DOUBLE PRECISION B( LDB, * ), D( * ), DL( * ), DU( * ) - function gtsv!(dl::StridedVector{$elty}, d::StridedVector{$elty}, du::StridedVector{$elty}, - B::StridedVecOrMat{$elty}) + function gtsv!(dl::AbstractVector{$elty}, d::AbstractVector{$elty}, du::AbstractVector{$elty}, + B::AbstractVecOrMat{$elty}) chkstride1(B, dl, d, du) n = length(d) if !(n >= length(dl) >= n - 1) @@ -2357,7 +2358,7 @@ for (gtsv, gttrf, gttrs, elty) in # .. Array Arguments .. # INTEGER IPIV( * ) # DOUBLE PRECISION D( * ), DL( * ), DU( * ), DU2( * ) - function gttrf!(dl::StridedVector{$elty}, d::StridedVector{$elty}, du::StridedVector{$elty}) + function gttrf!(dl::AbstractVector{$elty}, d::AbstractVector{$elty}, du::AbstractVector{$elty}) chkstride1(dl,d,du) n = length(d) if length(dl) != n - 1 @@ -2384,9 +2385,9 @@ for (gtsv, gttrf, gttrs, elty) in # .. Array Arguments .. # INTEGER IPIV( * ) # DOUBLE PRECISION B( LDB, * ), D( * ), DL( * ), DU( * ), DU2( * ) - function gttrs!(trans::Char, dl::StridedVector{$elty}, d::StridedVector{$elty}, - du::StridedVector{$elty}, du2::StridedVector{$elty}, ipiv::StridedVector{BlasInt}, - B::StridedVecOrMat{$elty}) + function gttrs!(trans::Char, dl::AbstractVector{$elty}, d::AbstractVector{$elty}, + du::AbstractVector{$elty}, du2::AbstractVector{$elty}, ipiv::AbstractVector{BlasInt}, + B::AbstractVecOrMat{$elty}) chktrans(trans) chkstride1(B, ipiv, dl, d, du, du2) n = length(d) @@ -2420,7 +2421,7 @@ superdiagonal. Overwrites `B` with the solution `X` and returns it. """ -gtsv!(dl::StridedVector, d::StridedVector, du::StridedVector, B::StridedVecOrMat) +gtsv!(dl::AbstractVector, d::AbstractVector, du::AbstractVector, B::AbstractVecOrMat) """ gttrf!(dl, d, du) -> (dl, d, du, du2, ipiv) @@ -2431,7 +2432,7 @@ subdiagonal, `d` on the diagonal, and `du` on the superdiagonal. Modifies `dl`, `d`, and `du` in-place and returns them and the second superdiagonal `du2` and the pivoting vector `ipiv`. """ -gttrf!(dl::StridedVector, d::StridedVector, du::StridedVector) +gttrf!(dl::AbstractVector, d::AbstractVector, du::AbstractVector) """ gttrs!(trans, dl, d, du, du2, ipiv, B) @@ -2440,8 +2441,8 @@ Solves the equation `A * X = B` (`trans = N`), `Transpose(A) * X = B` (`trans = or `Adjoint(A) * X = B` (`trans = C`) using the `LU` factorization computed by `gttrf!`. `B` is overwritten with the solution `X`. """ -gttrs!(trans::Char, dl::StridedVector, d::StridedVector, du::StridedVector, du2::StridedVector, - ipiv::StridedVector{BlasInt}, B::StridedVecOrMat) +gttrs!(trans::Char, dl::AbstractVector, d::AbstractVector, du::AbstractVector, du2::AbstractVector, + ipiv::AbstractVector{BlasInt}, B::AbstractVecOrMat) ## (OR) orthogonal (or UN, unitary) matrices, extractors and multiplication for (orglq, orgqr, orgql, orgrq, ormlq, ormqr, ormql, ormrq, gemqrt, elty) in @@ -2455,7 +2456,7 @@ for (orglq, orgqr, orgql, orgrq, ormlq, ormqr, ormql, ormrq, gemqrt, elty) in # INTEGER INFO, K, LDA, LWORK, M, N # * .. Array Arguments .. # DOUBLE PRECISION A( LDA, * ), TAU( * ), WORK( * ) - function orglq!(A::StridedMatrix{$elty}, tau::StridedVector{$elty}, k::Integer = length(tau)) + function orglq!(A::AbstractMatrix{$elty}, tau::AbstractVector{$elty}, k::Integer = length(tau)) chkstride1(A,tau) n = size(A, 2) m = min(n, size(A, 1)) @@ -2488,7 +2489,7 @@ for (orglq, orgqr, orgql, orgrq, ormlq, ormqr, ormql, ormrq, gemqrt, elty) in # INTEGER INFO, K, LDA, LWORK, M, N # * .. Array Arguments .. # DOUBLE PRECISION A( LDA, * ), TAU( * ), WORK( * ) - function orgqr!(A::StridedMatrix{$elty}, tau::StridedVector{$elty}, k::Integer = length(tau)) + function orgqr!(A::AbstractMatrix{$elty}, tau::AbstractVector{$elty}, k::Integer = length(tau)) chkstride1(A,tau) m = size(A, 1) n = min(m, size(A, 2)) @@ -2523,7 +2524,7 @@ for (orglq, orgqr, orgql, orgrq, ormlq, ormqr, ormql, ormrq, gemqrt, elty) in # INTEGER INFO, K, LDA, LWORK, M, N # * .. Array Arguments .. # DOUBLE PRECISION A( LDA, * ), TAU( * ), WORK( * ) - function orgql!(A::StridedMatrix{$elty}, tau::StridedVector{$elty}, k::Integer = length(tau)) + function orgql!(A::AbstractMatrix{$elty}, tau::AbstractVector{$elty}, k::Integer = length(tau)) chkstride1(A,tau) m = size(A, 1) n = min(m, size(A, 2)) @@ -2558,7 +2559,7 @@ for (orglq, orgqr, orgql, orgrq, ormlq, ormqr, ormql, ormrq, gemqrt, elty) in # INTEGER INFO, K, LDA, LWORK, M, N # * .. Array Arguments .. # DOUBLE PRECISION A( LDA, * ), TAU( * ), WORK( * ) - function orgrq!(A::StridedMatrix{$elty}, tau::StridedVector{$elty}, k::Integer = length(tau)) + function orgrq!(A::AbstractMatrix{$elty}, tau::AbstractVector{$elty}, k::Integer = length(tau)) chkstride1(A,tau) m, n = size(A) if n < m @@ -2593,8 +2594,8 @@ for (orglq, orgqr, orgql, orgrq, ormlq, ormqr, ormql, ormrq, gemqrt, elty) in # INTEGER INFO, K, LDA, LDC, LWORK, M, N # .. Array Arguments .. # DOUBLE PRECISION A( LDA, * ), C( LDC, * ), TAU( * ), WORK( * ) - function ormlq!(side::Char, trans::Char, A::StridedMatrix{$elty}, - tau::StridedVector{$elty}, C::StridedVecOrMat{$elty}) + function ormlq!(side::Char, trans::Char, A::AbstractMatrix{$elty}, + tau::AbstractVector{$elty}, C::AbstractVecOrMat{$elty}) chktrans(trans) chkside(side) chkstride1(A, C, tau) @@ -2639,8 +2640,8 @@ for (orglq, orgqr, orgql, orgrq, ormlq, ormqr, ormql, ormrq, gemqrt, elty) in # INTEGER INFO, K, LDA, LDC, M, N # .. Array Arguments .. # DOUBLE PRECISION A( LDA, * ), C( LDC, * ), TAU( * ), WORK( * ) - function ormqr!(side::Char, trans::Char, A::StridedMatrix{$elty}, - tau::StridedVector{$elty}, C::StridedVecOrMat{$elty}) + function ormqr!(side::Char, trans::Char, A::AbstractMatrix{$elty}, + tau::AbstractVector{$elty}, C::AbstractVecOrMat{$elty}) chktrans(trans) chkside(side) chkstride1(A, C, tau) @@ -2688,8 +2689,8 @@ for (orglq, orgqr, orgql, orgrq, ormlq, ormqr, ormql, ormrq, gemqrt, elty) in # INTEGER INFO, K, LDA, LDC, M, N # .. Array Arguments .. # DOUBLE PRECISION A( LDA, * ), C( LDC, * ), TAU( * ), WORK( * ) - function ormql!(side::Char, trans::Char, A::StridedMatrix{$elty}, - tau::StridedVector{$elty}, C::StridedVecOrMat{$elty}) + function ormql!(side::Char, trans::Char, A::AbstractMatrix{$elty}, + tau::AbstractVector{$elty}, C::AbstractVecOrMat{$elty}) chktrans(trans) chkside(side) chkstride1(A, C, tau) @@ -2737,8 +2738,8 @@ for (orglq, orgqr, orgql, orgrq, ormlq, ormqr, ormql, ormrq, gemqrt, elty) in # INTEGER INFO, K, LDA, LDC, LWORK, M, N # .. Array Arguments .. # DOUBLE PRECISION A( LDA, * ), C( LDC, * ), TAU( * ), WORK( * ) - function ormrq!(side::Char, trans::Char, A::StridedMatrix{$elty}, - tau::StridedVector{$elty}, C::StridedVecOrMat{$elty}) + function ormrq!(side::Char, trans::Char, A::AbstractMatrix{$elty}, + tau::AbstractVector{$elty}, C::AbstractVecOrMat{$elty}) chktrans(trans) chkside(side) chkstride1(A, C, tau) @@ -2776,7 +2777,7 @@ for (orglq, orgqr, orgql, orgrq, ormlq, ormqr, ormql, ormrq, gemqrt, elty) in C end - function gemqrt!(side::Char, trans::Char, V::StridedMatrix{$elty}, T::StridedMatrix{$elty}, C::StridedVecOrMat{$elty}) + function gemqrt!(side::Char, trans::Char, V::AbstractMatrix{$elty}, T::AbstractMatrix{$elty}, C::AbstractVecOrMat{$elty}) chktrans(trans) chkside(side) chkstride1(V, T, C) @@ -2837,7 +2838,7 @@ end Explicitly finds the matrix `Q` of a `LQ` factorization after calling `gelqf!` on `A`. Uses the output of `gelqf!`. `A` is overwritten by `Q`. """ -orglq!(A::StridedMatrix, tau::StridedVector, k::Integer = length(tau)) +orglq!(A::AbstractMatrix, tau::AbstractVector, k::Integer = length(tau)) """ orgqr!(A, tau, k = length(tau)) @@ -2845,7 +2846,7 @@ orglq!(A::StridedMatrix, tau::StridedVector, k::Integer = length(tau)) Explicitly finds the matrix `Q` of a `QR` factorization after calling `geqrf!` on `A`. Uses the output of `geqrf!`. `A` is overwritten by `Q`. """ -orgqr!(A::StridedMatrix, tau::StridedVector, k::Integer = length(tau)) +orgqr!(A::AbstractMatrix, tau::AbstractVector, k::Integer = length(tau)) """ orgql!(A, tau, k = length(tau)) @@ -2853,7 +2854,7 @@ orgqr!(A::StridedMatrix, tau::StridedVector, k::Integer = length(tau)) Explicitly finds the matrix `Q` of a `QL` factorization after calling `geqlf!` on `A`. Uses the output of `geqlf!`. `A` is overwritten by `Q`. """ -orgql!(A::StridedMatrix, tau::StridedVector, k::Integer = length(tau)) +orgql!(A::AbstractMatrix, tau::AbstractVector, k::Integer = length(tau)) """ orgrq!(A, tau, k = length(tau)) @@ -2861,7 +2862,7 @@ orgql!(A::StridedMatrix, tau::StridedVector, k::Integer = length(tau)) Explicitly finds the matrix `Q` of a `RQ` factorization after calling `gerqf!` on `A`. Uses the output of `gerqf!`. `A` is overwritten by `Q`. """ -orgrq!(A::StridedMatrix, tau::StridedVector, k::Integer = length(tau)) +orgrq!(A::AbstractMatrix, tau::AbstractVector, k::Integer = length(tau)) """ ormlq!(side, trans, A, tau, C) @@ -2871,7 +2872,7 @@ Computes `Q * C` (`trans = N`), `Transpose(Q) * C` (`trans = T`), `Adjoint(Q) * for `side = R` using `Q` from a `LQ` factorization of `A` computed using `gelqf!`. `C` is overwritten. """ -ormlq!(side::Char, trans::Char, A::StridedMatrix, tau::StridedVector, C::StridedVecOrMat) +ormlq!(side::Char, trans::Char, A::AbstractMatrix, tau::AbstractVector, C::AbstractVecOrMat) """ ormqr!(side, trans, A, tau, C) @@ -2881,7 +2882,7 @@ Computes `Q * C` (`trans = N`), `Transpose(Q) * C` (`trans = T`), `Adjoint(Q) * for `side = R` using `Q` from a `QR` factorization of `A` computed using `geqrf!`. `C` is overwritten. """ -ormqr!(side::Char, trans::Char, A::StridedMatrix, tau::StridedVector, C::StridedVecOrMat) +ormqr!(side::Char, trans::Char, A::AbstractMatrix, tau::AbstractVector, C::AbstractVecOrMat) """ ormql!(side, trans, A, tau, C) @@ -2891,7 +2892,7 @@ Computes `Q * C` (`trans = N`), `Transpose(Q) * C` (`trans = T`), `Adjoint(Q) * for `side = R` using `Q` from a `QL` factorization of `A` computed using `geqlf!`. `C` is overwritten. """ -ormql!(side::Char, trans::Char, A::StridedMatrix, tau::StridedVector, C::StridedVecOrMat) +ormql!(side::Char, trans::Char, A::AbstractMatrix, tau::AbstractVector, C::AbstractVecOrMat) """ ormrq!(side, trans, A, tau, C) @@ -2901,7 +2902,7 @@ Computes `Q * C` (`trans = N`), `Transpose(Q) * C` (`trans = T`), `Adjoint(Q) * for `side = R` using `Q` from a `RQ` factorization of `A` computed using `gerqf!`. `C` is overwritten. """ -ormrq!(side::Char, trans::Char, A::StridedMatrix, tau::StridedVector, C::StridedVecOrMat) +ormrq!(side::Char, trans::Char, A::AbstractMatrix, tau::AbstractVector, C::AbstractVecOrMat) """ gemqrt!(side, trans, V, T, C) @@ -2911,7 +2912,7 @@ Computes `Q * C` (`trans = N`), `Transpose(Q) * C` (`trans = T`), `Adjoint(Q) * for `side = R` using `Q` from a `QR` factorization of `A` computed using `geqrt!`. `C` is overwritten. """ -gemqrt!(side::Char, trans::Char, V::StridedMatrix, T::StridedMatrix, C::StridedVecOrMat) +gemqrt!(side::Char, trans::Char, V::AbstractMatrix, T::AbstractMatrix, C::AbstractVecOrMat) # (PO) positive-definite symmetric matrices, for (posv, potrf, potri, potrs, pstrf, elty, rtyp) in @@ -2926,7 +2927,7 @@ for (posv, potrf, potri, potrs, pstrf, elty, rtyp) in # INTEGER INFO, LDA, LDB, N, NRHS # .. Array Arguments .. # DOUBLE PRECISION A( LDA, * ), B( LDB, * ) - function posv!(uplo::Char, A::StridedMatrix{$elty}, B::StridedVecOrMat{$elty}) + function posv!(uplo::Char, A::AbstractMatrix{$elty}, B::AbstractVecOrMat{$elty}) chkstride1(A, B) n = checksquare(A) chkuplo(uplo) @@ -2949,7 +2950,7 @@ for (posv, potrf, potri, potrs, pstrf, elty, rtyp) in # INTEGER INFO, LDA, N # * .. Array Arguments .. # DOUBLE PRECISION A( LDA, * ) - function potrf!(uplo::Char, A::StridedMatrix{$elty}) + function potrf!(uplo::Char, A::AbstractMatrix{$elty}) chkstride1(A) checksquare(A) chkuplo(uplo) @@ -2974,7 +2975,7 @@ for (posv, potrf, potri, potrs, pstrf, elty, rtyp) in # INTEGER INFO, LDA, N # .. Array Arguments .. # DOUBLE PRECISION A( LDA, * ) - function potri!(uplo::Char, A::StridedMatrix{$elty}) + function potri!(uplo::Char, A::AbstractMatrix{$elty}) chkstride1(A) chkuplo(uplo) info = Ref{BlasInt}() @@ -2992,7 +2993,7 @@ for (posv, potrf, potri, potrs, pstrf, elty, rtyp) in # INTEGER INFO, LDA, LDB, N, NRHS # .. Array Arguments .. # DOUBLE PRECISION A( LDA, * ), B( LDB, * ) - function potrs!(uplo::Char, A::StridedMatrix{$elty}, B::StridedVecOrMat{$elty}) + function potrs!(uplo::Char, A::AbstractMatrix{$elty}, B::AbstractVecOrMat{$elty}) chkstride1(A, B) n = checksquare(A) chkuplo(uplo) @@ -3023,7 +3024,7 @@ for (posv, potrf, potri, potrs, pstrf, elty, rtyp) in # .. Array Arguments .. # DOUBLE PRECISION A( LDA, * ), WORK( 2*N ) # INTEGER PIV( N ) - function pstrf!(uplo::Char, A::StridedMatrix{$elty}, tol::Real) + function pstrf!(uplo::Char, A::AbstractMatrix{$elty}, tol::Real) chkstride1(A) n = checksquare(A) chkuplo(uplo) @@ -3050,7 +3051,7 @@ of `A` is computed. If `uplo = L` the lower Cholesky decomposition of `A` is computed. `A` is overwritten by its Cholesky decomposition. `B` is overwritten with the solution `X`. """ -posv!(uplo::Char, A::StridedMatrix, B::StridedVecOrMat) +posv!(uplo::Char, A::AbstractMatrix, B::AbstractVecOrMat) """ potrf!(uplo, A) @@ -3059,7 +3060,7 @@ Computes the Cholesky (upper if `uplo = U`, lower if `uplo = L`) decomposition of positive-definite matrix `A`. `A` is overwritten and returned with an info code. """ -potrf!(uplo::Char, A::StridedMatrix) +potrf!(uplo::Char, A::AbstractMatrix) """ potri!(uplo, A) @@ -3070,7 +3071,7 @@ decomposition. `A` is overwritten by its inverse and returned. """ -potri!(uplo::Char, A::StridedMatrix) +potri!(uplo::Char, A::AbstractMatrix) """ potrs!(uplo, A, B) @@ -3081,7 +3082,7 @@ positive definite matrix whose Cholesky decomposition was computed by computed. If `uplo = L` the lower Cholesky decomposition of `A` was computed. `B` is overwritten with the solution `X`. """ -potrs!(uplo::Char, A::StridedMatrix, B::StridedVecOrMat) +potrs!(uplo::Char, A::AbstractMatrix, B::AbstractVecOrMat) """ pstrf!(uplo, A, tol) -> (A, piv, rank, info) @@ -3094,7 +3095,7 @@ Returns `A`, the pivots `piv`, the rank of `A`, and an `info` code. If `info = 0 the factorization succeeded. If `info = i > 0 `, then `A` is indefinite or rank-deficient. """ -pstrf!(uplo::Char, A::StridedMatrix, tol::Real) +pstrf!(uplo::Char, A::AbstractMatrix, tol::Real) # (PT) positive-definite, symmetric, tri-diagonal matrices # Direct solvers for general tridiagonal and symmetric positive-definite tridiagonal @@ -3109,7 +3110,7 @@ for (ptsv, pttrf, elty, relty) in # INTEGER INFO, LDB, N, NRHS # .. Array Arguments .. # DOUBLE PRECISION B( LDB, * ), D( * ), E( * ) - function ptsv!(D::StridedVector{$relty}, E::StridedVector{$elty}, B::StridedVecOrMat{$elty}) + function ptsv!(D::AbstractVector{$relty}, E::AbstractVector{$elty}, B::AbstractVecOrMat{$elty}) chkstride1(B, D, E) n = length(D) if length(E) != n - 1 @@ -3132,7 +3133,7 @@ for (ptsv, pttrf, elty, relty) in # INTEGER INFO, N # .. Array Arguments .. # DOUBLE PRECISION D( * ), E( * ) - function pttrf!(D::StridedVector{$relty}, E::StridedVector{$elty}) + function pttrf!(D::AbstractVector{$relty}, E::AbstractVector{$elty}) chkstride1(D, E) n = length(D) if length(E) != n - 1 @@ -3155,7 +3156,7 @@ Solves `A * X = B` for positive-definite tridiagonal `A`. `D` is the diagonal of `A` and `E` is the off-diagonal. `B` is overwritten with the solution `X` and returned. """ -ptsv!(D::StridedVector, E::StridedVector, B::StridedVecOrMat) +ptsv!(D::AbstractVector, E::AbstractVector, B::AbstractVecOrMat) """ pttrf!(D, E) @@ -3164,7 +3165,7 @@ Computes the LDLt factorization of a positive-definite tridiagonal matrix with `D` as diagonal and `E` as off-diagonal. `D` and `E` are overwritten and returned. """ -pttrf!(D::StridedVector, E::StridedVector) +pttrf!(D::AbstractVector, E::AbstractVector) for (pttrs, elty, relty) in ((:dpttrs_,:Float64,:Float64), @@ -3175,7 +3176,7 @@ for (pttrs, elty, relty) in # INTEGER INFO, LDB, N, NRHS # .. Array Arguments .. # DOUBLE PRECISION B( LDB, * ), D( * ), E( * ) - function pttrs!(D::StridedVector{$relty}, E::StridedVector{$elty}, B::StridedVecOrMat{$elty}) + function pttrs!(D::AbstractVector{$relty}, E::AbstractVector{$elty}, B::AbstractVecOrMat{$elty}) chkstride1(B, D, E) n = length(D) if length(E) != n - 1 @@ -3207,7 +3208,7 @@ for (pttrs, elty, relty) in # * .. Array Arguments .. # DOUBLE PRECISION D( * ) # COMPLEX*16 B( LDB, * ), E( * ) - function pttrs!(uplo::Char, D::StridedVector{$relty}, E::StridedVector{$elty}, B::StridedVecOrMat{$elty}) + function pttrs!(uplo::Char, D::AbstractVector{$relty}, E::AbstractVector{$elty}, B::AbstractVecOrMat{$elty}) chkstride1(B, D, E) chkuplo(uplo) n = length(D) @@ -3235,7 +3236,7 @@ Solves `A * X = B` for positive-definite tridiagonal `A` with diagonal `D` and off-diagonal `E` after computing `A`'s LDLt factorization using `pttrf!`. `B` is overwritten with the solution `X`. """ -pttrs!(D::StridedVector, E::StridedVector, B::StridedVecOrMat) +pttrs!(D::AbstractVector, E::AbstractVector, B::AbstractVecOrMat) ## (TR) triangular matrices: solver and inverse for (trtri, trtrs, elty) in @@ -3250,7 +3251,7 @@ for (trtri, trtrs, elty) in # INTEGER INFO, LDA, N # .. Array Arguments .. # DOUBLE PRECISION A( LDA, * ) - function trtri!(uplo::Char, diag::Char, A::StridedMatrix{$elty}) + function trtri!(uplo::Char, diag::Char, A::AbstractMatrix{$elty}) chkstride1(A) n = checksquare(A) chkuplo(uplo) @@ -3272,7 +3273,7 @@ for (trtri, trtrs, elty) in # * .. Array Arguments .. # DOUBLE PRECISION A( LDA, * ), B( LDB, * ) function trtrs!(uplo::Char, trans::Char, diag::Char, - A::StridedMatrix{$elty}, B::StridedVecOrMat{$elty}) + A::AbstractMatrix{$elty}, B::AbstractVecOrMat{$elty}) chktrans(trans) chkdiag(diag) chkstride1(A) @@ -3301,7 +3302,7 @@ triangular matrix `A`. If `diag = N`, `A` has non-unit diagonal elements. If `diag = U`, all diagonal elements of `A` are one. `A` is overwritten with its inverse. """ -trtri!(uplo::Char, diag::Char, A::StridedMatrix) +trtri!(uplo::Char, diag::Char, A::AbstractMatrix) """ trtrs!(uplo, trans, diag, A, B) @@ -3312,7 +3313,7 @@ triangular matrix `A`. If `diag = N`, `A` has non-unit diagonal elements. If `diag = U`, all diagonal elements of `A` are one. `B` is overwritten with the solution `X`. """ -trtrs!(uplo::Char, trans::Char, diag::Char, A::StridedMatrix, B::StridedVecOrMat) +trtrs!(uplo::Char, trans::Char, diag::Char, A::AbstractMatrix, B::AbstractVecOrMat) #Eigenvector computation and condition number estimation for (trcon, trevc, trrfs, elty) in @@ -3328,7 +3329,7 @@ for (trcon, trevc, trrfs, elty) in # .. Array Arguments .. # INTEGER IWORK( * ) # DOUBLE PRECISION A( LDA, * ), WORK( * ) - function trcon!(norm::Char, uplo::Char, diag::Char, A::StridedMatrix{$elty}) + function trcon!(norm::Char, uplo::Char, diag::Char, A::AbstractMatrix{$elty}) chkstride1(A) chkdiag(diag) n = checksquare(A) @@ -3357,9 +3358,9 @@ for (trcon, trevc, trrfs, elty) in # LOGICAL SELECT( * ) # DOUBLE PRECISION T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ), #$ WORK( * ) - function trevc!(side::Char, howmny::Char, select::StridedVector{BlasInt}, T::StridedMatrix{$elty}, - VL::StridedMatrix{$elty} = similar(T), - VR::StridedMatrix{$elty} = similar(T)) + function trevc!(side::Char, howmny::Char, select::AbstractVector{BlasInt}, T::AbstractMatrix{$elty}, + VL::AbstractMatrix{$elty} = similar(T), + VR::AbstractMatrix{$elty} = similar(T)) # Extract if side ∉ ['L','R','B'] throw(ArgumentError("side argument must be 'L' (left eigenvectors), 'R' (right eigenvectors), or 'B' (both), got $side")) @@ -3368,7 +3369,7 @@ for (trcon, trevc, trrfs, elty) in ldt, ldvl, ldvr = stride(T, 2), stride(VL, 2), stride(VR, 2) # Check - chkstride1(T, select) + chkstride1(T, select, VL, VR) # Allocate m = Ref{BlasInt}() @@ -3416,9 +3417,10 @@ for (trcon, trevc, trrfs, elty) in # DOUBLE PRECISION A( LDA, * ), B( LDB, * ), BERR( * ), FERR( * ), #$ WORK( * ), X( LDX, * ) function trrfs!(uplo::Char, trans::Char, diag::Char, - A::StridedMatrix{$elty}, B::StridedVecOrMat{$elty}, X::StridedVecOrMat{$elty}, - Ferr::StridedVector{$elty} = similar(B, $elty, size(B,2)), - Berr::StridedVector{$elty} = similar(B, $elty, size(B,2))) + A::AbstractMatrix{$elty}, B::AbstractVecOrMat{$elty}, X::AbstractVecOrMat{$elty}, + Ferr::AbstractVector{$elty} = similar(B, $elty, size(B,2)), + Berr::AbstractVector{$elty} = similar(B, $elty, size(B,2))) + chkstride1(A, B, X, Ferr, Berr) chktrans(trans) chkuplo(uplo) chkdiag(diag) @@ -3456,7 +3458,7 @@ for (trcon, trevc, trrfs, elty, relty) in # .. Array Arguments .. # DOUBLE PRECISION RWORK( * ) # COMPLEX*16 A( LDA, * ), WORK( * ) - function trcon!(norm::Char, uplo::Char, diag::Char, A::StridedMatrix{$elty}) + function trcon!(norm::Char, uplo::Char, diag::Char, A::AbstractMatrix{$elty}) chkstride1(A) n = checksquare(A) chkuplo(uplo) @@ -3486,15 +3488,15 @@ for (trcon, trevc, trrfs, elty, relty) in # DOUBLE PRECISION RWORK( * ) # COMPLEX*16 T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ), #$ WORK( * ) - function trevc!(side::Char, howmny::Char, select::StridedVector{BlasInt}, T::StridedMatrix{$elty}, - VL::StridedMatrix{$elty} = similar(T), - VR::StridedMatrix{$elty} = similar(T)) + function trevc!(side::Char, howmny::Char, select::AbstractVector{BlasInt}, T::AbstractMatrix{$elty}, + VL::AbstractMatrix{$elty} = similar(T), + VR::AbstractMatrix{$elty} = similar(T)) # Extract n, mm = checksquare(T), size(VL, 2) ldt, ldvl, ldvr = stride(T, 2), stride(VL, 2), stride(VR, 2) # Check - chkstride1(T, select) + chkstride1(T, select, VL, VR) if side ∉ ['L','R','B'] throw(ArgumentError("side argument must be 'L' (left eigenvectors), 'R' (right eigenvectors), or 'B' (both), got $side")) end @@ -3545,9 +3547,10 @@ for (trcon, trevc, trrfs, elty, relty) in # DOUBLE PRECISION A( LDA, * ), B( LDB, * ), BERR( * ), FERR( * ), #$ WORK( * ), X( LDX, * ) function trrfs!(uplo::Char, trans::Char, diag::Char, - A::StridedMatrix{$elty}, B::StridedVecOrMat{$elty}, X::StridedVecOrMat{$elty}, - Ferr::StridedVector{$relty} = similar(B, $relty, size(B,2)), - Berr::StridedVector{$relty} = similar(B, $relty, size(B,2))) + A::AbstractMatrix{$elty}, B::AbstractVecOrMat{$elty}, X::AbstractVecOrMat{$elty}, + Ferr::AbstractVector{$relty} = similar(B, $relty, size(B,2)), + Berr::AbstractVector{$relty} = similar(B, $relty, size(B,2))) + chkstride1(A, B, X, Ferr, Berr) chktrans(trans) chkuplo(uplo) chkdiag(diag) @@ -3581,7 +3584,7 @@ diagonal elements. If `diag = U`, all diagonal elements of `A` are one. If `norm = I`, the condition number is found in the infinity norm. If `norm = O` or `1`, the condition number is found in the one norm. """ -trcon!(norm::Char, uplo::Char, diag::Char, A::StridedMatrix) +trcon!(norm::Char, uplo::Char, diag::Char, A::AbstractMatrix) """ trevc!(side, howmny, select, T, VL = similar(T), VR = similar(T)) @@ -3594,8 +3597,8 @@ eigenvectors are found and backtransformed using `VL` and `VR`. If `howmny = S`, only the eigenvectors corresponding to the values in `select` are computed. """ -trevc!(side::Char, howmny::Char, select::StridedVector{BlasInt}, T::StridedMatrix, - VL::StridedMatrix = similar(T), VR::StridedMatrix = similar(T)) +trevc!(side::Char, howmny::Char, select::AbstractVector{BlasInt}, T::AbstractMatrix, + VL::AbstractMatrix = similar(T), VR::AbstractMatrix = similar(T)) """ trrfs!(uplo, trans, diag, A, B, X, Ferr, Berr) -> (Ferr, Berr) @@ -3609,8 +3612,8 @@ diagonal elements. If `diag = U`, all diagonal elements of `A` are one. `Ferr` and `Berr` are optional inputs. `Ferr` is the forward error and `Berr` is the backward error, each component-wise. """ -trrfs!(uplo::Char, trans::Char, diag::Char, A::StridedMatrix, B::StridedVecOrMat, - X::StridedVecOrMat, Ferr::StridedVector, Berr::StridedVector) +trrfs!(uplo::Char, trans::Char, diag::Char, A::AbstractMatrix, B::AbstractVecOrMat, + X::AbstractVecOrMat, Ferr::AbstractVector, Berr::AbstractVector) ## (ST) Symmetric tridiagonal - eigendecomposition for (stev, stebz, stegr, stein, elty) in @@ -3620,7 +3623,7 @@ for (stev, stebz, stegr, stein, elty) in # , (:cstev_,:ComplexF32) ) @eval begin - function stev!(job::Char, dv::StridedVector{$elty}, ev::StridedVector{$elty}) + function stev!(job::Char, dv::AbstractVector{$elty}, ev::AbstractVector{$elty}) chkstride1(dv, ev) n = length(dv) if length(ev) != n - 1 @@ -3641,7 +3644,7 @@ for (stev, stebz, stegr, stein, elty) in #* matrix T. The user may ask for all eigenvalues, all eigenvalues #* in the half-open interval (VL, VU], or the IL-th through IU-th #* eigenvalues. - function stebz!(range::Char, order::Char, vl::$elty, vu::$elty, il::Integer, iu::Integer, abstol::Real, dv::StridedVector{$elty}, ev::StridedVector{$elty}) + function stebz!(range::Char, order::Char, vl::$elty, vu::$elty, il::Integer, iu::Integer, abstol::Real, dv::AbstractVector{$elty}, ev::AbstractVector{$elty}) chkstride1(dv, ev) n = length(dv) if length(ev) != n - 1 @@ -3671,7 +3674,7 @@ for (stev, stebz, stegr, stein, elty) in w[1:m[]], iblock[1:m[]], isplit[1:nsplit[1]] end - function stegr!(jobz::Char, range::Char, dv::StridedVector{$elty}, ev::StridedVector{$elty}, vl::Real, vu::Real, il::Integer, iu::Integer) + function stegr!(jobz::Char, range::Char, dv::AbstractVector{$elty}, ev::AbstractVector{$elty}, vl::Real, vu::Real, il::Integer, iu::Integer) chkstride1(dv, ev) n = length(dv) if length(ev) != n - 1 @@ -3712,7 +3715,7 @@ for (stev, stebz, stegr, stein, elty) in m[] == length(w) ? w : w[1:m[]], m[] == size(Z, 2) ? Z : Z[:,1:m[]] end - function stein!(dv::StridedVector{$elty}, ev_in::StridedVector{$elty}, w_in::StridedVector{$elty}, iblock_in::StridedVector{BlasInt}, isplit_in::StridedVector{BlasInt}) + function stein!(dv::AbstractVector{$elty}, ev_in::AbstractVector{$elty}, w_in::AbstractVector{$elty}, iblock_in::AbstractVector{BlasInt}, isplit_in::AbstractVector{BlasInt}) chkstride1(dv, ev_in, w_in, iblock_in, isplit_in) n = length(dv) if length(ev_in) != n - 1 @@ -3762,12 +3765,12 @@ for (stev, stebz, stegr, stein, elty) in end end end -stegr!(jobz::Char, dv::StridedVector, ev::StridedVector) = stegr!(jobz, 'A', dv, ev, 0.0, 0.0, 0, 0) +stegr!(jobz::Char, dv::AbstractVector, ev::AbstractVector) = stegr!(jobz, 'A', dv, ev, 0.0, 0.0, 0, 0) # Allow user to skip specification of iblock and isplit -stein!(dv::StridedVector, ev::StridedVector, w_in::StridedVector)=stein!(dv, ev, w_in, zeros(BlasInt,0), zeros(BlasInt,0)) +stein!(dv::AbstractVector, ev::AbstractVector, w_in::AbstractVector) = stein!(dv, ev, w_in, zeros(BlasInt,0), zeros(BlasInt,0)) # Allow user to specify just one eigenvector to get in stein! -stein!(dv::StridedVector, ev::StridedVector, eval::Real)=stein!(dv, ev, [eval], zeros(BlasInt,0), zeros(BlasInt,0)) +stein!(dv::AbstractVector, ev::AbstractVector, eval::Real) = stein!(dv, ev, [eval], zeros(BlasInt,0), zeros(BlasInt,0)) """ stev!(job, dv, ev) -> (dv, Zmat) @@ -3777,7 +3780,7 @@ diagonal and `ev` as off-diagonal. If `job = N` only the eigenvalues are found and returned in `dv`. If `job = V` then the eigenvectors are also found and returned in `Zmat`. """ -stev!(job::Char, dv::StridedVector, ev::StridedVector) +stev!(job::Char, dv::AbstractVector, ev::AbstractVector) """ stebz!(range, order, vl, vu, il, iu, abstol, dv, ev) -> (dv, iblock, isplit) @@ -3790,7 +3793,7 @@ are found. If `range = V`, the eigenvalues in the half-open interval block. If `order = E`, they are ordered across all the blocks. `abstol` can be set as a tolerance for convergence. """ -stebz!(range::Char, order::Char, vl, vu, il::Integer, iu::Integer, abstol::Real, dv::StridedVector, ev::StridedVector) +stebz!(range::Char, order::Char, vl, vu, il::Integer, iu::Integer, abstol::Real, dv::AbstractVector, ev::AbstractVector) """ stegr!(jobz, range, dv, ev, vl, vu, il, iu) -> (w, Z) @@ -3803,7 +3806,7 @@ are found. If `range = V`, the eigenvalues in the half-open interval `il` and `iu` are found. The eigenvalues are returned in `w` and the eigenvectors in `Z`. """ -stegr!(jobz::Char, range::Char, dv::StridedVector, ev::StridedVector, vl::Real, vu::Real, il::Integer, iu::Integer) +stegr!(jobz::Char, range::Char, dv::AbstractVector, ev::AbstractVector, vl::Real, vu::Real, il::Integer, iu::Integer) """ stein!(dv, ev_in, w_in, iblock_in, isplit_in) @@ -3814,7 +3817,7 @@ eigenvalues for which to find corresponding eigenvectors. `iblock_in` specifies the submatrices corresponding to the eigenvalues in `w_in`. `isplit_in` specifies the splitting points between the submatrix blocks. """ -stein!(dv::StridedVector, ev_in::StridedVector, w_in::StridedVector, iblock_in::StridedVector{BlasInt}, isplit_in::StridedVector{BlasInt}) +stein!(dv::AbstractVector, ev_in::AbstractVector, w_in::AbstractVector, iblock_in::AbstractVector{BlasInt}, isplit_in::AbstractVector{BlasInt}) ## (SY) symmetric real matrices - Bunch-Kaufman decomposition, ## solvers (direct and factored) and inverse. @@ -3829,7 +3832,7 @@ for (syconv, sysv, sytrf, sytri, sytrs, elty) in # * .. Array Arguments .. # INTEGER IPIV( * ) # DOUBLE PRECISION A( LDA, * ), WORK( * ) - function syconv!(uplo::Char, A::StridedMatrix{$elty}, ipiv::StridedVector{BlasInt}) + function syconv!(uplo::Char, A::AbstractMatrix{$elty}, ipiv::AbstractVector{BlasInt}) chkstride1(A, ipiv) n = checksquare(A) chkuplo(uplo) @@ -3851,7 +3854,7 @@ for (syconv, sysv, sytrf, sytri, sytrs, elty) in # .. Array Arguments .. # INTEGER IPIV( * ) # DOUBLE PRECISION A( LDA, * ), B( LDB, * ), WORK( * ) - function sysv!(uplo::Char, A::StridedMatrix{$elty}, B::StridedVecOrMat{$elty}) + function sysv!(uplo::Char, A::AbstractMatrix{$elty}, B::AbstractVecOrMat{$elty}) chkstride1(A,B) n = checksquare(A) chkuplo(uplo) @@ -3885,7 +3888,7 @@ for (syconv, sysv, sytrf, sytri, sytrs, elty) in # * .. Array Arguments .. # INTEGER IPIV( * ) # DOUBLE PRECISION A( LDA, * ), WORK( * ) - function sytrf!(uplo::Char, A::StridedMatrix{$elty}) + function sytrf!(uplo::Char, A::AbstractMatrix{$elty}) chkstride1(A) n = checksquare(A) chkuplo(uplo) @@ -3917,7 +3920,7 @@ for (syconv, sysv, sytrf, sytri, sytrs, elty) in # * .. Array Arguments .. # INTEGER IPIV( * ) # DOUBLE PRECISION A( LDA, * ), WORK( * ) -# function sytri!(uplo::Char, A::StridedMatrix{$elty}, ipiv::Vector{BlasInt}) +# function sytri!(uplo::Char, A::AbstractMatrix{$elty}, ipiv::Vector{BlasInt}) # chkstride1(A) # n = checksquare(A) # chkuplo(uplo) @@ -3946,7 +3949,7 @@ for (syconv, sysv, sytrf, sytri, sytrs, elty) in # .. Array Arguments .. # INTEGER IPIV( * ) # DOUBLE PRECISION A( LDA, * ), WORK( * ) - function sytri!(uplo::Char, A::StridedMatrix{$elty}, ipiv::StridedVector{BlasInt}) + function sytri!(uplo::Char, A::AbstractMatrix{$elty}, ipiv::AbstractVector{BlasInt}) chkstride1(A, ipiv) n = checksquare(A) chkuplo(uplo) @@ -3969,8 +3972,8 @@ for (syconv, sysv, sytrf, sytri, sytrs, elty) in # .. Array Arguments .. # INTEGER IPIV( * ) # DOUBLE PRECISION A( LDA, * ), B( LDB, * ) - function sytrs!(uplo::Char, A::StridedMatrix{$elty}, - ipiv::StridedVector{BlasInt}, B::StridedVecOrMat{$elty}) + function sytrs!(uplo::Char, A::AbstractMatrix{$elty}, + ipiv::AbstractVector{BlasInt}, B::AbstractVecOrMat{$elty}) chkstride1(A,B,ipiv) n = checksquare(A) chkuplo(uplo) @@ -4001,7 +4004,7 @@ for (sysv, sytrf, sytri, sytrs, syconvf, elty) in # .. Array Arguments .. # INTEGER IPIV( * ) # DOUBLE PRECISION A( LDA, * ), B( LDB, * ), WORK( * ) - function sysv_rook!(uplo::Char, A::StridedMatrix{$elty}, B::StridedVecOrMat{$elty}) + function sysv_rook!(uplo::Char, A::AbstractMatrix{$elty}, B::AbstractVecOrMat{$elty}) chkstride1(A,B) n = checksquare(A) chkuplo(uplo) @@ -4035,7 +4038,7 @@ for (sysv, sytrf, sytri, sytrs, syconvf, elty) in # * .. Array Arguments .. # INTEGER IPIV( * ) # DOUBLE PRECISION A( LDA, * ), WORK( * ) - function sytrf_rook!(uplo::Char, A::StridedMatrix{$elty}) + function sytrf_rook!(uplo::Char, A::AbstractMatrix{$elty}) chkstride1(A) n = checksquare(A) chkuplo(uplo) @@ -4067,7 +4070,7 @@ for (sysv, sytrf, sytri, sytrs, syconvf, elty) in # .. Array Arguments .. # INTEGER IPIV( * ) # DOUBLE PRECISION A( LDA, * ), WORK( * ) - function sytri_rook!(uplo::Char, A::StridedMatrix{$elty}, ipiv::StridedVector{BlasInt}) + function sytri_rook!(uplo::Char, A::AbstractMatrix{$elty}, ipiv::AbstractVector{BlasInt}) chkstride1(A, ipiv) n = checksquare(A) chkuplo(uplo) @@ -4090,8 +4093,8 @@ for (sysv, sytrf, sytri, sytrs, syconvf, elty) in # .. Array Arguments .. # INTEGER IPIV( * ) # DOUBLE PRECISION A( LDA, * ), B( LDB, * ) - function sytrs_rook!(uplo::Char, A::StridedMatrix{$elty}, - ipiv::StridedVector{BlasInt}, B::StridedVecOrMat{$elty}) + function sytrs_rook!(uplo::Char, A::AbstractMatrix{$elty}, + ipiv::AbstractVector{BlasInt}, B::AbstractVecOrMat{$elty}) chkstride1(A,B,ipiv) n = checksquare(A) chkuplo(uplo) @@ -4117,8 +4120,8 @@ for (sysv, sytrf, sytri, sytrs, syconvf, elty) in # INTEGER IPIV( * ) # DOUBLE PRECISION A( LDA, * ), E( * ) function syconvf_rook!(uplo::Char, way::Char, - A::StridedMatrix{$elty}, ipiv::StridedVector{BlasInt}, - e::StridedVector{$elty} = Vector{$elty}(uninitialized, length(ipiv))) + A::AbstractMatrix{$elty}, ipiv::AbstractVector{BlasInt}, + e::AbstractVector{$elty} = Vector{$elty}(uninitialized, length(ipiv))) # extract n = checksquare(A) lda = max(1, stride(A, 2)) @@ -4165,7 +4168,7 @@ for (syconv, hesv, hetrf, hetri, hetrs, elty, relty) in # .. Array Arguments .. # INTEGER IPIV( * ) # COMPLEX*16 A( LDA, * ), WORK( * ) - function syconv!(uplo::Char, A::StridedMatrix{$elty}, ipiv::StridedVector{BlasInt}) + function syconv!(uplo::Char, A::AbstractMatrix{$elty}, ipiv::AbstractVector{BlasInt}) chkstride1(A,ipiv) n = checksquare(A) chkuplo(uplo) @@ -4187,7 +4190,7 @@ for (syconv, hesv, hetrf, hetri, hetrs, elty, relty) in # * .. Array Arguments .. # INTEGER IPIV( * ) # COMPLEX*16 A( LDA, * ), B( LDB, * ), WORK( * ) - function hesv!(uplo::Char, A::StridedMatrix{$elty}, B::StridedVecOrMat{$elty}) + function hesv!(uplo::Char, A::AbstractMatrix{$elty}, B::AbstractVecOrMat{$elty}) chkstride1(A,B) n = checksquare(A) chkuplo(uplo) @@ -4221,7 +4224,7 @@ for (syconv, hesv, hetrf, hetri, hetrs, elty, relty) in # * .. Array Arguments .. # INTEGER IPIV( * ) # COMPLEX*16 A( LDA, * ), WORK( * ) - function hetrf!(uplo::Char, A::StridedMatrix{$elty}) + function hetrf!(uplo::Char, A::AbstractMatrix{$elty}) chkstride1(A) n = checksquare(A) chkuplo(uplo) @@ -4251,7 +4254,7 @@ for (syconv, hesv, hetrf, hetri, hetrs, elty, relty) in # * .. Array Arguments .. # INTEGER IPIV( * ) # COMPLEX*16 A( LDA, * ), WORK( * ) -# function hetri!(uplo::Char, A::StridedMatrix{$elty}, ipiv::Vector{BlasInt}) +# function hetri!(uplo::Char, A::AbstractMatrix{$elty}, ipiv::Vector{BlasInt}) # chkstride1(A) # n = checksquare(A) # chkuplo(uplo) @@ -4281,7 +4284,7 @@ for (syconv, hesv, hetrf, hetri, hetrs, elty, relty) in # * .. Array Arguments .. # INTEGER IPIV( * ) # COMPLEX*16 A( LDA, * ), WORK( * ) - function hetri!(uplo::Char, A::StridedMatrix{$elty}, ipiv::StridedVector{BlasInt}) + function hetri!(uplo::Char, A::AbstractMatrix{$elty}, ipiv::AbstractVector{BlasInt}) chkstride1(A, ipiv) n = checksquare(A) chkuplo(uplo) @@ -4303,8 +4306,8 @@ for (syconv, hesv, hetrf, hetri, hetrs, elty, relty) in # * .. Array Arguments .. # INTEGER IPIV( * ) # COMPLEX*16 A( LDA, * ), B( LDB, * ) - function hetrs!(uplo::Char, A::StridedMatrix{$elty}, - ipiv::StridedVector{BlasInt}, B::StridedVecOrMat{$elty}) + function hetrs!(uplo::Char, A::AbstractMatrix{$elty}, + ipiv::AbstractVector{BlasInt}, B::AbstractVecOrMat{$elty}) chkstride1(A,B,ipiv) n = checksquare(A) if n != size(B,1) @@ -4333,7 +4336,7 @@ for (hesv, hetrf, hetri, hetrs, elty, relty) in # * .. Array Arguments .. # INTEGER IPIV( * ) # COMPLEX*16 A( LDA, * ), B( LDB, * ), WORK( * ) - function hesv_rook!(uplo::Char, A::StridedMatrix{$elty}, B::StridedVecOrMat{$elty}) + function hesv_rook!(uplo::Char, A::AbstractMatrix{$elty}, B::AbstractVecOrMat{$elty}) chkstride1(A,B) n = checksquare(A) chkuplo(uplo) @@ -4367,7 +4370,7 @@ for (hesv, hetrf, hetri, hetrs, elty, relty) in # * .. Array Arguments .. # INTEGER IPIV( * ) # COMPLEX*16 A( LDA, * ), WORK( * ) - function hetrf_rook!(uplo::Char, A::StridedMatrix{$elty}) + function hetrf_rook!(uplo::Char, A::AbstractMatrix{$elty}) chkstride1(A) n = checksquare(A) chkuplo(uplo) @@ -4397,7 +4400,7 @@ for (hesv, hetrf, hetri, hetrs, elty, relty) in # * .. Array Arguments .. # INTEGER IPIV( * ) # COMPLEX*16 A( LDA, * ), WORK( * ) - function hetri_rook!(uplo::Char, A::StridedMatrix{$elty}, ipiv::StridedVector{BlasInt}) + function hetri_rook!(uplo::Char, A::AbstractMatrix{$elty}, ipiv::AbstractVector{BlasInt}) chkstride1(A,ipiv) n = checksquare(A) chkuplo(uplo) @@ -4419,8 +4422,8 @@ for (hesv, hetrf, hetri, hetrs, elty, relty) in # * .. Array Arguments .. # INTEGER IPIV( * ) # COMPLEX*16 A( LDA, * ), B( LDB, * ) - function hetrs_rook!(uplo::Char, A::StridedMatrix{$elty}, - ipiv::StridedVector{BlasInt}, B::StridedVecOrMat{$elty}) + function hetrs_rook!(uplo::Char, A::AbstractMatrix{$elty}, + ipiv::AbstractVector{BlasInt}, B::AbstractVecOrMat{$elty}) chkstride1(A,B,ipiv) n = checksquare(A) if n != size(B,1) @@ -4450,7 +4453,7 @@ for (sysv, sytrf, sytri, sytrs, elty, relty) in # * .. Array Arguments .. # INTEGER IPIV( * ) # COMPLEX*16 A( LDA, * ), B( LDB, * ), WORK( * ) - function sysv!(uplo::Char, A::StridedMatrix{$elty}, B::StridedVecOrMat{$elty}) + function sysv!(uplo::Char, A::AbstractMatrix{$elty}, B::AbstractVecOrMat{$elty}) chkstride1(A,B) n = checksquare(A) chkuplo(uplo) @@ -4485,7 +4488,7 @@ for (sysv, sytrf, sytri, sytrs, elty, relty) in # * .. Array Arguments .. # INTEGER IPIV( * ) # COMPLEX*16 A( LDA, * ), WORK( * ) - function sytrf!(uplo::Char, A::StridedMatrix{$elty}) + function sytrf!(uplo::Char, A::AbstractMatrix{$elty}) chkstride1(A) n = checksquare(A) chkuplo(uplo) @@ -4518,7 +4521,7 @@ for (sysv, sytrf, sytri, sytrs, elty, relty) in # * .. Array Arguments .. # INTEGER IPIV( * ) # COMPLEX*16 A( LDA, * ), WORK( * ) -# function sytri!(uplo::Char, A::StridedMatrix{$elty}, ipiv::Vector{BlasInt}) +# function sytri!(uplo::Char, A::AbstractMatrix{$elty}, ipiv::Vector{BlasInt}) # chkstride1(A) # n = checksquare(A) # chkuplo(uplo) @@ -4547,7 +4550,7 @@ for (sysv, sytrf, sytri, sytrs, elty, relty) in # * .. Array Arguments .. # INTEGER IPIV( * ) # COMPLEX*16 A( LDA, * ), WORK( * ) - function sytri!(uplo::Char, A::StridedMatrix{$elty}, ipiv::StridedVector{BlasInt}) + function sytri!(uplo::Char, A::AbstractMatrix{$elty}, ipiv::AbstractVector{BlasInt}) chkstride1(A, ipiv) n = checksquare(A) chkuplo(uplo) @@ -4569,8 +4572,8 @@ for (sysv, sytrf, sytri, sytrs, elty, relty) in # * .. Array Arguments .. # INTEGER IPIV( * ) # COMPLEX*16 A( LDA, * ), B( LDB, * ) - function sytrs!(uplo::Char, A::StridedMatrix{$elty}, - ipiv::StridedVector{BlasInt}, B::StridedVecOrMat{$elty}) + function sytrs!(uplo::Char, A::AbstractMatrix{$elty}, + ipiv::AbstractVector{BlasInt}, B::AbstractVecOrMat{$elty}) chkstride1(A,B,ipiv) n = checksquare(A) chkuplo(uplo) @@ -4601,7 +4604,7 @@ for (sysv, sytrf, sytri, sytrs, syconvf, elty, relty) in # * .. Array Arguments .. # INTEGER IPIV( * ) # COMPLEX*16 A( LDA, * ), B( LDB, * ), WORK( * ) - function sysv_rook!(uplo::Char, A::StridedMatrix{$elty}, B::StridedVecOrMat{$elty}) + function sysv_rook!(uplo::Char, A::AbstractMatrix{$elty}, B::AbstractVecOrMat{$elty}) chkstride1(A,B) n = checksquare(A) chkuplo(uplo) @@ -4636,7 +4639,7 @@ for (sysv, sytrf, sytri, sytrs, syconvf, elty, relty) in # * .. Array Arguments .. # INTEGER IPIV( * ) # COMPLEX*16 A( LDA, * ), WORK( * ) - function sytrf_rook!(uplo::Char, A::StridedMatrix{$elty}) + function sytrf_rook!(uplo::Char, A::AbstractMatrix{$elty}) chkstride1(A) n = checksquare(A) chkuplo(uplo) @@ -4669,7 +4672,7 @@ for (sysv, sytrf, sytri, sytrs, syconvf, elty, relty) in # * .. Array Arguments .. # INTEGER IPIV( * ) # COMPLEX*16 A( LDA, * ), WORK( * ) - function sytri_rook!(uplo::Char, A::StridedMatrix{$elty}, ipiv::StridedVector{BlasInt}) + function sytri_rook!(uplo::Char, A::AbstractMatrix{$elty}, ipiv::AbstractVector{BlasInt}) chkstride1(A, ipiv) n = checksquare(A) chkuplo(uplo) @@ -4691,8 +4694,8 @@ for (sysv, sytrf, sytri, sytrs, syconvf, elty, relty) in # * .. Array Arguments .. # INTEGER IPIV( * ) # COMPLEX*16 A( LDA, * ), B( LDB, * ) - function sytrs_rook!(uplo::Char, A::StridedMatrix{$elty}, - ipiv::StridedVector{BlasInt}, B::StridedVecOrMat{$elty}) + function sytrs_rook!(uplo::Char, A::AbstractMatrix{$elty}, + ipiv::AbstractVector{BlasInt}, B::AbstractVecOrMat{$elty}) chkstride1(A,B,ipiv) n = checksquare(A) chkuplo(uplo) @@ -4718,8 +4721,10 @@ for (sysv, sytrf, sytri, sytrs, syconvf, elty, relty) in # INTEGER IPIV( * ) # COMPLEX*16 A( LDA, * ), E( * ) function syconvf_rook!(uplo::Char, way::Char, - A::StridedMatrix{$elty}, ipiv::StridedVector{BlasInt}, - e::StridedVector{$elty} = Vector{$elty}(uninitialized, length(ipiv))) + A::AbstractMatrix{$elty}, ipiv::AbstractVector{BlasInt}, + e::AbstractVector{$elty} = Vector{$elty}(uninitialized, length(ipiv))) + chkstride1(A, ipiv, e) + # extract n = checksquare(A) lda = stride(A, 2) @@ -4760,7 +4765,7 @@ is upper triangular. If `uplo = L`, it is lower triangular. `ipiv` is the pivot vector from the triangular factorization. `A` is overwritten by `L` and `D`. """ -syconv!(uplo::Char, A::StridedMatrix, ipiv::StridedVector{BlasInt}) +syconv!(uplo::Char, A::AbstractMatrix, ipiv::AbstractVector{BlasInt}) """ sysv!(uplo, A, B) -> (B, A, ipiv) @@ -4771,7 +4776,7 @@ the upper half of `A` is stored. If `uplo = L`, the lower half is stored. Bunch-Kaufman factorization. `ipiv` contains pivoting information about the factorization. """ -sysv!(uplo::Char, A::StridedMatrix, B::StridedVecOrMat) +sysv!(uplo::Char, A::AbstractMatrix, B::AbstractVecOrMat) """ sytrf!(uplo, A) -> (A, ipiv, info) @@ -4785,7 +4790,7 @@ the error code `info` which is a non-negative integer. If `info` is positive the matrix is singular and the diagonal part of the factorization is exactly zero at position `info`. """ -sytrf!(uplo::Char, A::StridedMatrix) +sytrf!(uplo::Char, A::AbstractMatrix) """ sytri!(uplo, A, ipiv) @@ -4794,7 +4799,7 @@ Computes the inverse of a symmetric matrix `A` using the results of `sytrf!`. If `uplo = U`, the upper half of `A` is stored. If `uplo = L`, the lower half is stored. `A` is overwritten by its inverse. """ -sytri!(uplo::Char, A::StridedMatrix, ipiv::StridedVector{BlasInt}) +sytri!(uplo::Char, A::AbstractMatrix, ipiv::AbstractVector{BlasInt}) """ sytrs!(uplo, A, ipiv, B) @@ -4804,7 +4809,7 @@ results of `sytrf!`. If `uplo = U`, the upper half of `A` is stored. If `uplo = L`, the lower half is stored. `B` is overwritten by the solution `X`. """ -sytrs!(uplo::Char, A::StridedMatrix, ipiv::StridedVector{BlasInt}, B::StridedVecOrMat) +sytrs!(uplo::Char, A::AbstractMatrix, ipiv::AbstractVector{BlasInt}, B::AbstractVecOrMat) """ @@ -4816,7 +4821,7 @@ the upper half of `A` is stored. If `uplo = L`, the lower half is stored. Bunch-Kaufman factorization. `ipiv` contains pivoting information about the factorization. """ -hesv!(uplo::Char, A::StridedMatrix, B::StridedVecOrMat) +hesv!(uplo::Char, A::AbstractMatrix, B::AbstractVecOrMat) """ hetrf!(uplo, A) -> (A, ipiv, info) @@ -4830,7 +4835,7 @@ the error code `info` which is a non-negative integer. If `info` is positive the matrix is singular and the diagonal part of the factorization is exactly zero at position `info`. """ -hetrf!(uplo::Char, A::StridedMatrix) +hetrf!(uplo::Char, A::AbstractMatrix) """ hetri!(uplo, A, ipiv) @@ -4839,7 +4844,7 @@ Computes the inverse of a Hermitian matrix `A` using the results of `sytrf!`. If `uplo = U`, the upper half of `A` is stored. If `uplo = L`, the lower half is stored. `A` is overwritten by its inverse. """ -hetri!(uplo::Char, A::StridedMatrix, ipiv::StridedVector{BlasInt}) +hetri!(uplo::Char, A::AbstractMatrix, ipiv::AbstractVector{BlasInt}) """ hetrs!(uplo, A, ipiv, B) @@ -4849,7 +4854,7 @@ results of `sytrf!`. If `uplo = U`, the upper half of `A` is stored. If `uplo = L`, the lower half is stored. `B` is overwritten by the solution `X`. """ -hetrs!(uplo::Char, A::StridedMatrix, ipiv::StridedVector{BlasInt}, B::StridedVecOrMat) +hetrs!(uplo::Char, A::AbstractMatrix, ipiv::AbstractVector{BlasInt}, B::AbstractVecOrMat) # Symmetric (real) eigensolvers for (syev, syevr, sygvd, elty) in @@ -4862,7 +4867,7 @@ for (syev, syevr, sygvd, elty) in # INTEGER INFO, LDA, LWORK, N # * .. Array Arguments .. # DOUBLE PRECISION A( LDA, * ), W( * ), WORK( * ) - function syev!(jobz::Char, uplo::Char, A::StridedMatrix{$elty}) + function syev!(jobz::Char, uplo::Char, A::AbstractMatrix{$elty}) chkstride1(A) n = checksquare(A) W = similar(A, $elty, n) @@ -4894,7 +4899,7 @@ for (syev, syevr, sygvd, elty) in # * .. Array Arguments .. # INTEGER ISUPPZ( * ), IWORK( * ) # DOUBLE PRECISION A( LDA, * ), W( * ), WORK( * ), Z( LDZ, * ) - function syevr!(jobz::Char, range::Char, uplo::Char, A::StridedMatrix{$elty}, + function syevr!(jobz::Char, range::Char, uplo::Char, A::AbstractMatrix{$elty}, vl::AbstractFloat, vu::AbstractFloat, il::Integer, iu::Integer, abstol::AbstractFloat) chkstride1(A) n = checksquare(A) @@ -4943,7 +4948,7 @@ for (syev, syevr, sygvd, elty) in end w[1:m[]], Z[:,1:(jobz == 'V' ? m[] : 0)] end - syevr!(jobz::Char, A::StridedMatrix{$elty}) = + syevr!(jobz::Char, A::AbstractMatrix{$elty}) = syevr!(jobz, 'A', 'U', A, 0.0, 0.0, 0, 0, -1.0) # Generalized eigenproblem @@ -4956,7 +4961,7 @@ for (syev, syevr, sygvd, elty) in # * .. Array Arguments .. # INTEGER IWORK( * ) # DOUBLE PRECISION A( LDA, * ), B( LDB, * ), W( * ), WORK( * ) - function sygvd!(itype::Integer, jobz::Char, uplo::Char, A::StridedMatrix{$elty}, B::StridedMatrix{$elty}) + function sygvd!(itype::Integer, jobz::Char, uplo::Char, A::AbstractMatrix{$elty}, B::AbstractMatrix{$elty}) chkstride1(A, B) n, m = checksquare(A, B) if n != m @@ -5006,7 +5011,7 @@ for (syev, syevr, sygvd, elty, relty) in # * .. Array Arguments .. # DOUBLE PRECISION RWORK( * ), W( * ) # COMPLEX*16 A( LDA, * ), WORK( * ) - function syev!(jobz::Char, uplo::Char, A::StridedMatrix{$elty}) + function syev!(jobz::Char, uplo::Char, A::AbstractMatrix{$elty}) chkstride1(A) n = checksquare(A) W = similar(A, $relty, n) @@ -5041,7 +5046,7 @@ for (syev, syevr, sygvd, elty, relty) in # INTEGER ISUPPZ( * ), IWORK( * ) # DOUBLE PRECISION RWORK( * ), W( * ) # COMPLEX*16 A( LDA, * ), WORK( * ), Z( LDZ, * ) - function syevr!(jobz::Char, range::Char, uplo::Char, A::StridedMatrix{$elty}, + function syevr!(jobz::Char, range::Char, uplo::Char, A::AbstractMatrix{$elty}, vl::AbstractFloat, vu::AbstractFloat, il::Integer, iu::Integer, abstol::AbstractFloat) chkstride1(A) n = checksquare(A) @@ -5095,7 +5100,7 @@ for (syev, syevr, sygvd, elty, relty) in end w[1:m[]], Z[:,1:(jobz == 'V' ? m[] : 0)] end - syevr!(jobz::Char, A::StridedMatrix{$elty}) = + syevr!(jobz::Char, A::AbstractMatrix{$elty}) = syevr!(jobz, 'A', 'U', A, 0.0, 0.0, 0, 0, -1.0) # SUBROUTINE ZHEGVD( ITYPE, JOBZ, UPLO, N, A, LDA, B, LDB, W, WORK, @@ -5108,7 +5113,7 @@ for (syev, syevr, sygvd, elty, relty) in # INTEGER IWORK( * ) # DOUBLE PRECISION RWORK( * ), W( * ) # COMPLEX*16 A( LDA, * ), B( LDB, * ), WORK( * ) - function sygvd!(itype::Integer, jobz::Char, uplo::Char, A::StridedMatrix{$elty}, B::StridedMatrix{$elty}) + function sygvd!(itype::Integer, jobz::Char, uplo::Char, A::AbstractMatrix{$elty}, B::AbstractMatrix{$elty}) chkstride1(A, B) n, m = checksquare(A, B) if n != m @@ -5157,7 +5162,7 @@ Finds the eigenvalues (`jobz = N`) or eigenvalues and eigenvectors (`jobz = V`) of a symmetric matrix `A`. If `uplo = U`, the upper triangle of `A` is used. If `uplo = L`, the lower triangle of `A` is used. """ -syev!(jobz::Char, uplo::Char, A::StridedMatrix) +syev!(jobz::Char, uplo::Char, A::AbstractMatrix) """ syevr!(jobz, range, uplo, A, vl, vu, il, iu, abstol) -> (W, Z) @@ -5172,7 +5177,7 @@ found. `abstol` can be set as a tolerance for convergence. The eigenvalues are returned in `W` and the eigenvectors in `Z`. """ -syevr!(jobz::Char, range::Char, uplo::Char, A::StridedMatrix, +syevr!(jobz::Char, range::Char, uplo::Char, A::AbstractMatrix, vl::AbstractFloat, vu::AbstractFloat, il::Integer, iu::Integer, abstol::AbstractFloat) """ @@ -5187,7 +5192,7 @@ of `A` and `B` are used. If `uplo = L`, the lower triangles of `A` and `A * B * x = lambda * x`. If `itype = 3`, the problem to solve is `B * A * x = lambda * x`. """ -sygvd!(itype::Integer, jobz::Char, uplo::Char, A::StridedMatrix, B::StridedMatrix) +sygvd!(itype::Integer, jobz::Char, uplo::Char, A::AbstractMatrix, B::AbstractMatrix) ## (BD) Bidiagonal matrices - singular value decomposition for (bdsqr, relty, elty) in @@ -5196,9 +5201,9 @@ for (bdsqr, relty, elty) in (:zbdsqr_,:Float64,:ComplexF64), (:cbdsqr_,:Float32,:ComplexF32)) @eval begin - function bdsqr!(uplo::Char, d::StridedVector{$relty}, e_::StridedVector{$relty}, - Vt::StridedMatrix{$elty}, U::StridedMatrix{$elty}, C::StridedMatrix{$elty}) - chkstride1(d, e_) + function bdsqr!(uplo::Char, d::AbstractVector{$relty}, e_::AbstractVector{$relty}, + Vt::AbstractMatrix{$elty}, U::AbstractMatrix{$elty}, C::AbstractMatrix{$elty}) + chkstride1(d, e_, Vt, U, C) # Extract number n = length(d) ncvt, nru, ncc = size(Vt, 2), size(U, 1), size(C, 2) @@ -5248,7 +5253,7 @@ compute the product `Q' * C`. Returns the singular values in `d`, and the matrix `C` overwritten with `Q' * C`. """ -bdsqr!(uplo::Char, d::StridedVector, e_::StridedVector, Vt::StridedMatrix, U::StridedMatrix, C::StridedMatrix) +bdsqr!(uplo::Char, d::AbstractVector, e_::AbstractVector, Vt::AbstractMatrix, U::AbstractMatrix, C::AbstractMatrix) #Defined only for real types for (bdsdc, elty) in @@ -5266,7 +5271,7 @@ for (bdsdc, elty) in # INTEGER IQ( * ), IWORK( * ) # DOUBLE PRECISION D( * ), E( * ), Q( * ), U( LDU, * ), # $ VT( LDVT, * ), WORK( * ) - function bdsdc!(uplo::Char, compq::Char, d::StridedVector{$elty}, e_::StridedVector{$elty}) + function bdsdc!(uplo::Char, compq::Char, d::AbstractVector{$elty}, e_::AbstractVector{$elty}) chkstride1(d, e_) n, ldiq, ldq, ldu, ldvt = length(d), 1, 1, 1, 1 chkuplo(uplo) @@ -5319,7 +5324,7 @@ and vectors are found in compact form. Only works for real types. Returns the singular values in `d`, and if `compq = P`, the compact singular vectors in `iq`. """ -bdsdc!(uplo::Char, compq::Char, d::StridedVector, e_::StridedVector) +bdsdc!(uplo::Char, compq::Char, d::AbstractVector, e_::AbstractVector) for (gecon, elty) in ((:dgecon_,:Float64), @@ -5335,7 +5340,7 @@ for (gecon, elty) in # * .. Array Arguments .. # INTEGER IWORK( * ) # DOUBLE PRECISION A( LDA, * ), WORK( * ) - function gecon!(normtype::Char, A::StridedMatrix{$elty}, anorm::$elty) + function gecon!(normtype::Char, A::AbstractMatrix{$elty}, anorm::$elty) chkstride1(A) n = checksquare(A) lda = max(1, stride(A, 2)) @@ -5369,7 +5374,7 @@ for (gecon, elty, relty) in # * .. Array Arguments .. # DOUBLE PRECISION RWORK( * ) # COMPLEX*16 A( LDA, * ), WORK( * ) - function gecon!(normtype::Char, A::StridedMatrix{$elty}, anorm::$relty) + function gecon!(normtype::Char, A::AbstractMatrix{$elty}, anorm::$relty) chkstride1(A) n = checksquare(A) lda = max(1, stride(A, 2)) @@ -5397,7 +5402,7 @@ the condition number is found in the infinity norm. If `normtype = O` or `1`, the condition number is found in the one norm. `A` must be the result of `getrf!` and `anorm` is the norm of `A` in the relevant norm. """ -gecon!(normtype::Char, A::StridedMatrix, anorm) +gecon!(normtype::Char, A::AbstractMatrix, anorm) for (gehrd, elty) in ((:dgehrd_,:Float64), @@ -5411,7 +5416,7 @@ for (gehrd, elty) in # * .. # * .. Array Arguments .. # DOUBLE PRECISION A( LDA, * ), TAU( * ), WORK( * ) - function gehrd!(ilo::Integer, ihi::Integer, A::StridedMatrix{$elty}) + function gehrd!(ilo::Integer, ihi::Integer, A::AbstractMatrix{$elty}) chkstride1(A) n = checksquare(A) chkfinite(A) # balancing routines don't support NaNs and Infs @@ -5437,7 +5442,7 @@ for (gehrd, elty) in end end end -gehrd!(A::StridedMatrix) = gehrd!(1, size(A, 1), A) +gehrd!(A::AbstractMatrix) = gehrd!(1, size(A, 1), A) """ gehrd!(ilo, ihi, A) -> (A, tau) @@ -5447,7 +5452,7 @@ then `ilo` and `ihi` are the outputs of `gebal!`. Otherwise they should be `ilo = 1` and `ihi = size(A,2)`. `tau` contains the elementary reflectors of the factorization. """ -gehrd!(ilo::Integer, ihi::Integer, A::StridedMatrix) +gehrd!(ilo::Integer, ihi::Integer, A::AbstractMatrix) for (orghr, elty) in ((:dorghr_,:Float64), @@ -5460,7 +5465,7 @@ for (orghr, elty) in # * .. # * .. Array Arguments .. # DOUBLE PRECISION A( LDA, * ), TAU( * ), WORK( * ) - function orghr!(ilo::Integer, ihi::Integer, A::StridedMatrix{$elty}, tau::StridedVector{$elty}) + function orghr!(ilo::Integer, ihi::Integer, A::AbstractMatrix{$elty}, tau::AbstractVector{$elty}) chkstride1(A, tau) n = checksquare(A) if n - length(tau) != 1 @@ -5494,7 +5499,7 @@ end Explicitly finds `Q`, the orthogonal/unitary matrix from `gehrd!`. `ilo`, `ihi`, `A`, and `tau` must correspond to the input/output to `gehrd!`. """ -orghr!(ilo::Integer, ihi::Integer, A::StridedMatrix, tau::StridedVector) +orghr!(ilo::Integer, ihi::Integer, A::AbstractMatrix, tau::AbstractVector) for (ormhr, elty) in ((:dormhr_,:Float64), @@ -5508,10 +5513,10 @@ for (ormhr, elty) in # .. # .. Array Arguments .. # DOUBLE PRECISION a( lda, * ), c( ldc, * ), tau( * ), work( * ) - function ormhr!(side::Char, trans::Char, ilo::Integer, ihi::Integer, A::StridedMatrix{$elty}, - tau::StridedVector{$elty}, C::StridedVecOrMat{$elty}) + function ormhr!(side::Char, trans::Char, ilo::Integer, ihi::Integer, A::AbstractMatrix{$elty}, + tau::AbstractVector{$elty}, C::AbstractVecOrMat{$elty}) - chkstride1(A, tau) + chkstride1(A, tau, C) n = checksquare(A) mC, nC = size(C, 1), size(C, 2) @@ -5558,7 +5563,7 @@ for (gees, gges, elty) in # LOGICAL BWORK( * ) # DOUBLE PRECISION A( LDA, * ), VS( LDVS, * ), WI( * ), WORK( * ), # $ WR( * ) - function gees!(jobvs::Char, A::StridedMatrix{$elty}) + function gees!(jobvs::Char, A::AbstractMatrix{$elty}) chkstride1(A) n = checksquare(A) sdim = Vector{BlasInt}(uninitialized, 1) @@ -5597,7 +5602,7 @@ for (gees, gges, elty) in # DOUBLE PRECISION A( LDA, * ), ALPHAI( * ), ALPHAR( * ), # $ B( LDB, * ), BETA( * ), VSL( LDVSL, * ), # $ VSR( LDVSR, * ), WORK( * ) - function gges!(jobvsl::Char, jobvsr::Char, A::StridedMatrix{$elty}, B::StridedMatrix{$elty}) + function gges!(jobvsl::Char, jobvsr::Char, A::AbstractMatrix{$elty}, B::AbstractMatrix{$elty}) chkstride1(A, B) n, m = checksquare(A, B) if n != m @@ -5651,7 +5656,7 @@ for (gees, gges, elty, relty) in # LOGICAL BWORK( * ) # DOUBLE PRECISION RWORK( * ) # COMPLEX*16 A( LDA, * ), VS( LDVS, * ), W( * ), WORK( * ) - function gees!(jobvs::Char, A::StridedMatrix{$elty}) + function gees!(jobvs::Char, A::AbstractMatrix{$elty}) chkstride1(A) n = checksquare(A) sort = 'N' @@ -5692,7 +5697,7 @@ for (gees, gges, elty, relty) in # COMPLEX*16 A( LDA, * ), ALPHA( * ), B( LDB, * ), # $ BETA( * ), VSL( LDVSL, * ), VSR( LDVSR, * ), # $ WORK( * ) - function gges!(jobvsl::Char, jobvsr::Char, A::StridedMatrix{$elty}, B::StridedMatrix{$elty}) + function gges!(jobvsl::Char, jobvsr::Char, A::AbstractMatrix{$elty}, B::AbstractMatrix{$elty}) chkstride1(A, B) n, m = checksquare(A, B) if n != m @@ -5743,7 +5748,7 @@ vectors (`jobvs = V`) of matrix `A`. `A` is overwritten by its Schur form. Returns `A`, `vs` containing the Schur vectors, and `w`, containing the eigenvalues. """ -gees!(jobvs::Char, A::StridedMatrix) +gees!(jobvs::Char, A::AbstractMatrix) """ @@ -5756,7 +5761,7 @@ vectors (`jobsvl = V`), or right Schur vectors (`jobvsr = V`) of `A` and The generalized eigenvalues are returned in `alpha` and `beta`. The left Schur vectors are returned in `vsl` and the right Schur vectors are returned in `vsr`. """ -gges!(jobvsl::Char, jobvsr::Char, A::StridedMatrix, B::StridedMatrix) +gges!(jobvsl::Char, jobvsr::Char, A::AbstractMatrix, B::AbstractMatrix) for (trexc, trsen, tgsen, elty) in ((:dtrexc_, :dtrsen_, :dtgsen_, :Float64), @@ -5768,7 +5773,7 @@ for (trexc, trsen, tgsen, elty) in # * .. # * .. Array Arguments .. # DOUBLE PRECISION Q( LDQ, * ), T( LDT, * ), WORK( * ) - function trexc!(compq::Char, ifst::BlasInt, ilst::BlasInt, T::StridedMatrix{$elty}, Q::StridedMatrix{$elty}) + function trexc!(compq::Char, ifst::BlasInt, ilst::BlasInt, T::AbstractMatrix{$elty}, Q::AbstractMatrix{$elty}) chkstride1(T, Q) n = checksquare(T) ldt = max(1, stride(T, 2)) @@ -5787,7 +5792,7 @@ for (trexc, trsen, tgsen, elty) in chklapackerror(info[]) T, Q end - trexc!(ifst::BlasInt, ilst::BlasInt, T::StridedMatrix{$elty}, Q::StridedMatrix{$elty}) = + trexc!(ifst::BlasInt, ilst::BlasInt, T::AbstractMatrix{$elty}, Q::AbstractMatrix{$elty}) = trexc!('V', ifst, ilst, T, Q) # * .. Scalar Arguments .. @@ -5799,8 +5804,8 @@ for (trexc, trsen, tgsen, elty) in # LOGICAL SELECT( * ) # INTEGER IWORK( * ) # DOUBLE PRECISION Q( LDQ, * ), T( LDT, * ), WI( * ), WORK( * ), WR( * ) - function trsen!(job::Char, compq::Char, select::StridedVector{BlasInt}, - T::StridedMatrix{$elty}, Q::StridedMatrix{$elty}) + function trsen!(job::Char, compq::Char, select::AbstractVector{BlasInt}, + T::AbstractMatrix{$elty}, Q::AbstractMatrix{$elty}) chkstride1(T, Q, select) n = checksquare(T) ldt = max(1, stride(T, 2)) @@ -5838,7 +5843,7 @@ for (trexc, trsen, tgsen, elty) in end T, Q, iszero(wi) ? wr : complex.(wr, wi), s[], sep[] end - trsen!(select::StridedVector{BlasInt}, T::StridedMatrix{$elty}, Q::StridedMatrix{$elty}) = + trsen!(select::AbstractVector{BlasInt}, T::AbstractMatrix{$elty}, Q::AbstractMatrix{$elty}) = trsen!('N', 'V', select, T, Q) # .. Scalar Arguments .. @@ -5854,8 +5859,8 @@ for (trexc, trsen, tgsen, elty) in # $ B( LDB, * ), BETA( * ), DIF( * ), Q( LDQ, * ), # $ WORK( * ), Z( LDZ, * ) # .. - function tgsen!(select::StridedVector{BlasInt}, S::StridedMatrix{$elty}, T::StridedMatrix{$elty}, - Q::StridedMatrix{$elty}, Z::StridedMatrix{$elty}) + function tgsen!(select::AbstractVector{BlasInt}, S::AbstractMatrix{$elty}, T::AbstractMatrix{$elty}, + Q::AbstractMatrix{$elty}, Z::AbstractMatrix{$elty}) chkstride1(select, S, T, Q, Z) n, nt, nq, nz = checksquare(S, T, Q, Z) if n != nt @@ -5920,7 +5925,7 @@ for (trexc, trsen, tgsen, elty, relty) in # .. # .. Array Arguments .. # DOUBLE PRECISION Q( LDQ, * ), T( LDT, * ), WORK( * ) - function trexc!(compq::Char, ifst::BlasInt, ilst::BlasInt, T::StridedMatrix{$elty}, Q::StridedMatrix{$elty}) + function trexc!(compq::Char, ifst::BlasInt, ilst::BlasInt, T::AbstractMatrix{$elty}, Q::AbstractMatrix{$elty}) chkstride1(T, Q) n = checksquare(T) ldt = max(1, stride(T, 2)) @@ -5938,7 +5943,7 @@ for (trexc, trsen, tgsen, elty, relty) in chklapackerror(info[]) T, Q end - trexc!(ifst::BlasInt, ilst::BlasInt, T::StridedMatrix{$elty}, Q::StridedMatrix{$elty}) = + trexc!(ifst::BlasInt, ilst::BlasInt, T::AbstractMatrix{$elty}, Q::AbstractMatrix{$elty}) = trexc!('V', ifst, ilst, T, Q) # .. Scalar Arguments .. @@ -5949,8 +5954,8 @@ for (trexc, trsen, tgsen, elty, relty) in # .. Array Arguments .. # LOGICAL SELECT( * ) # COMPLEX Q( LDQ, * ), T( LDT, * ), W( * ), WORK( * ) - function trsen!(job::Char, compq::Char, select::StridedVector{BlasInt}, - T::StridedMatrix{$elty}, Q::StridedMatrix{$elty}) + function trsen!(job::Char, compq::Char, select::AbstractVector{BlasInt}, + T::AbstractMatrix{$elty}, Q::AbstractMatrix{$elty}) chkstride1(select, T, Q) n = checksquare(T) ldt = max(1, stride(T, 2)) @@ -5983,7 +5988,7 @@ for (trexc, trsen, tgsen, elty, relty) in end T, Q, w, s[], sep[] end - trsen!(select::StridedVector{BlasInt}, T::StridedMatrix{$elty}, Q::StridedMatrix{$elty}) = + trsen!(select::AbstractVector{BlasInt}, T::AbstractMatrix{$elty}, Q::AbstractMatrix{$elty}) = trsen!('N', 'V', select, T, Q) # .. Scalar Arguments .. @@ -5999,8 +6004,8 @@ for (trexc, trsen, tgsen, elty, relty) in # COMPLEX*16 A( LDA, * ), ALPHA( * ), B( LDB, * ), # $ BETA( * ), Q( LDQ, * ), WORK( * ), Z( LDZ, * ) # .. - function tgsen!(select::StridedVector{BlasInt}, S::StridedMatrix{$elty}, T::StridedMatrix{$elty}, - Q::StridedMatrix{$elty}, Z::StridedMatrix{$elty}) + function tgsen!(select::AbstractVector{BlasInt}, S::AbstractMatrix{$elty}, T::AbstractMatrix{$elty}, + Q::AbstractMatrix{$elty}, Z::AbstractMatrix{$elty}) chkstride1(select, S, T, Q, Z) n, nt, nq, nz = checksquare(S, T, Q, Z) if n != nt @@ -6061,7 +6066,7 @@ Reorder the Schur factorization of a matrix. If `compq = V`, the Schur vectors `Q` are reordered. If `compq = N` they are not modified. `ifst` and `ilst` specify the reordering of the vectors. """ -trexc!(compq::Char, ifst::BlasInt, ilst::BlasInt, T::StridedMatrix, Q::StridedMatrix) +trexc!(compq::Char, ifst::BlasInt, ilst::BlasInt, T::AbstractMatrix, Q::AbstractMatrix) """ trsen!(compq, job, select, T, Q) -> (T, Q, w, s, sep) @@ -6079,7 +6084,7 @@ Returns `T`, `Q`, reordered eigenvalues in `w`, the condition number of the cluster of eigenvalues `s`, and the condition number of the invariant subspace `sep`. """ -trsen!(compq::Char, job::Char, select::StridedVector{BlasInt}, T::StridedMatrix, Q::StridedMatrix) +trsen!(compq::Char, job::Char, select::AbstractVector{BlasInt}, T::AbstractMatrix, Q::AbstractMatrix) """ tgsen!(select, S, T, Q, Z) -> (S, T, alpha, beta, Q, Z) @@ -6087,15 +6092,15 @@ trsen!(compq::Char, job::Char, select::StridedVector{BlasInt}, T::StridedMatrix, Reorders the vectors of a generalized Schur decomposition. `select` specifices the eigenvalues in each cluster. """ -tgsen!(select::StridedVector{BlasInt}, S::StridedMatrix, T::StridedMatrix, Q::StridedMatrix, Z::StridedMatrix) +tgsen!(select::AbstractVector{BlasInt}, S::AbstractMatrix, T::AbstractMatrix, Q::AbstractMatrix, Z::AbstractMatrix) for (fn, elty, relty) in ((:dtrsyl_, :Float64, :Float64), (:strsyl_, :Float32, :Float32), (:ztrsyl_, :ComplexF64, :Float64), (:ctrsyl_, :ComplexF32, :Float32)) @eval begin - function trsyl!(transa::Char, transb::Char, A::StridedMatrix{$elty}, - B::StridedMatrix{$elty}, C::StridedMatrix{$elty}, isgn::Int=1) + function trsyl!(transa::Char, transb::Char, A::AbstractMatrix{$elty}, + B::AbstractMatrix{$elty}, C::AbstractMatrix{$elty}, isgn::Int=1) chkstride1(A, B, C) m, n = checksquare(A, B) lda = max(1, stride(A, 2)) @@ -6132,6 +6137,6 @@ transposed. Similarly for `transb` and `B`. If `isgn = 1`, the equation Returns `X` (overwriting `C`) and `scale`. """ -trsyl!(transa::Char, transb::Char, A::StridedMatrix, B::StridedMatrix, C::StridedMatrix, isgn::Int=1) +trsyl!(transa::Char, transb::Char, A::AbstractMatrix, B::AbstractMatrix, C::AbstractMatrix, isgn::Int=1) end # module diff --git a/base/linalg/linalg.jl b/base/linalg/linalg.jl index 90a617ddfb24b..130e35c3118d0 100644 --- a/base/linalg/linalg.jl +++ b/base/linalg/linalg.jl @@ -13,8 +13,8 @@ import Base: USE_BLAS64, abs, acos, acosh, acot, acoth, acsc, acsch, adjoint, as cosh, cot, coth, csc, csch, eltype, exp, findmax, findmin, fill!, floor, getindex, hcat, getproperty, imag, inv, isapprox, isone, IndexStyle, kron, length, log, map, ndims, oneunit, parent, power_by_squaring, print_matrix, promote_rule, real, round, sec, sech, - setindex!, show, similar, sin, sincos, sinh, size, sqrt, tan, tanh, transpose, trunc, - typed_hcat, vec + setindex!, show, similar, sin, sincos, sinh, size, size_to_strides, sqrt, StridedReinterpretArray, + StridedReshapedArray, strides, stride, tan, tanh, transpose, trunc, typed_hcat, vec using Base: hvcat_fill, iszero, IndexLinear, _length, promote_op, promote_typeof, @propagate_inbounds, @pure, reduce, typed_vcat # We use `_length` because of non-1 indices; releases after julia 0.5 @@ -161,9 +161,40 @@ end # Check that stride of matrix/vector is 1 # Writing like this to avoid splatting penalty when called with multiple arguments, # see PR 16416 +""" + stride1(A) -> Int + +Return the distance between successive array elements +in dimension 1 in units of element size. + +# Examples +```jldoctest +julia> A = [1,2,3,4] +4-element Array{Int64,1}: + 1 + 2 + 3 + 4 + +julia> Base.LinAlg.stride1(A) +1 + +julia> B = view(A, 2:2:4) +2-element view(::Array{Int64,1}, 2:2:4) with eltype Int64: + 2 + 4 + +julia> Base.LinAlg.stride1(B) +2 +``` +""" +stride1(x) = stride(x,1) +stride1(x::Array) = 1 +stride1(x::DenseArray) = stride(x, 1)::Int + @inline chkstride1(A...) = _chkstride1(true, A...) @noinline _chkstride1(ok::Bool) = ok || error("matrix does not have contiguous columns") -@inline _chkstride1(ok::Bool, A, B...) = _chkstride1(ok & (stride(A, 1) == 1), B...) +@inline _chkstride1(ok::Bool, A, B...) = _chkstride1(ok & (stride1(A) == 1), B...) """ LinAlg.checksquare(A) diff --git a/base/sparse/sparsevector.jl b/base/sparse/sparsevector.jl index 4d86349377166..489e7e2c0c366 100644 --- a/base/sparse/sparsevector.jl +++ b/base/sparse/sparsevector.jl @@ -1356,8 +1356,8 @@ for (vop, fun, mode) in [(:_vadd, :+, 1), (:_vmul, :*, 0)] @eval begin $(vop)(x::AbstractSparseVector, y::AbstractSparseVector) = _binarymap($(fun), x, y, $mode) - $(vop)(x::StridedVector, y::AbstractSparseVector) = _binarymap($(fun), x, y, $mode) - $(vop)(x::AbstractSparseVector, y::StridedVector) = _binarymap($(fun), x, y, $mode) + $(vop)(x::AbstractVector, y::AbstractSparseVector) = _binarymap($(fun), x, y, $mode) + $(vop)(x::AbstractSparseVector, y::AbstractVector) = _binarymap($(fun), x, y, $mode) end end @@ -1370,13 +1370,13 @@ broadcast(::typeof(*), x::AbstractSparseVector{Bool}, y::BitVector) = _vmul(x, y for (op, vop) in [(:+, :_vadd), (:-, :_vsub), (:*, :_vmul)] op != :* && @eval begin $(op)(x::AbstractSparseVector, y::AbstractSparseVector) = $(vop)(x, y) - $(op)(x::StridedVector, y::AbstractSparseVector) = $(vop)(x, y) - $(op)(x::AbstractSparseVector, y::StridedVector) = $(vop)(x, y) + $(op)(x::AbstractVector, y::AbstractSparseVector) = $(vop)(x, y) + $(op)(x::AbstractSparseVector, y::AbstractVector) = $(vop)(x, y) end @eval begin broadcast(::typeof($op), x::AbstractSparseVector, y::AbstractSparseVector) = $(vop)(x, y) - broadcast(::typeof($op), x::StridedVector, y::AbstractSparseVector) = $(vop)(x, y) - broadcast(::typeof($op), x::AbstractSparseVector, y::StridedVector) = $(vop)(x, y) + broadcast(::typeof($op), x::AbstractVector, y::AbstractSparseVector) = $(vop)(x, y) + broadcast(::typeof($op), x::AbstractSparseVector, y::AbstractVector) = $(vop)(x, y) end end @@ -1384,17 +1384,17 @@ end broadcast(::typeof(min), x::SparseVector{<:Real}, y::SparseVector{<:Real}) = _binarymap(min, x, y, 2) broadcast(::typeof(min), x::AbstractSparseVector{<:Real}, y::AbstractSparseVector{<:Real}) = _binarymap(min, x, y, 2) -broadcast(::typeof(min), x::StridedVector{<:Real}, y::AbstractSparseVector{<:Real}) = _binarymap(min, x, y, 2) -broadcast(::typeof(min), x::AbstractSparseVector{<:Real}, y::StridedVector{<:Real}) = _binarymap(min, x, y, 2) +broadcast(::typeof(min), x::AbstractVector{<:Real}, y::AbstractSparseVector{<:Real}) = _binarymap(min, x, y, 2) +broadcast(::typeof(min), x::AbstractSparseVector{<:Real}, y::AbstractVector{<:Real}) = _binarymap(min, x, y, 2) broadcast(::typeof(max), x::SparseVector{<:Real}, y::SparseVector{<:Real}) = _binarymap(max, x, y, 2) broadcast(::typeof(max), x::AbstractSparseVector{<:Real}, y::AbstractSparseVector{<:Real}) = _binarymap(max, x, y, 2) -broadcast(::typeof(max), x::StridedVector{<:Real}, y::AbstractSparseVector{<:Real}) = _binarymap(max, x, y, 2) -broadcast(::typeof(max), x::AbstractSparseVector{<:Real}, y::StridedVector{<:Real}) = _binarymap(max, x, y, 2) +broadcast(::typeof(max), x::AbstractVector{<:Real}, y::AbstractSparseVector{<:Real}) = _binarymap(max, x, y, 2) +broadcast(::typeof(max), x::AbstractSparseVector{<:Real}, y::AbstractVector{<:Real}) = _binarymap(max, x, y, 2) complex(x::AbstractSparseVector{<:Real}, y::AbstractSparseVector{<:Real}) = _binarymap(complex, x, y, 1) -complex(x::StridedVector{<:Real}, y::AbstractSparseVector{<:Real}) = _binarymap(complex, x, y, 1) -complex(x::AbstractSparseVector{<:Real}, y::StridedVector{<:Real}) = _binarymap(complex, x, y, 1) +complex(x::AbstractVector{<:Real}, y::AbstractSparseVector{<:Real}) = _binarymap(complex, x, y, 1) +complex(x::AbstractSparseVector{<:Real}, y::AbstractVector{<:Real}) = _binarymap(complex, x, y, 1) ### Reduction @@ -1440,7 +1440,7 @@ adjoint(sv::SparseVector) = Adjoint(sv) # axpy -function LinAlg.axpy!(a::Number, x::SparseVectorUnion, y::StridedVector) +function LinAlg.axpy!(a::Number, x::SparseVectorUnion, y::AbstractVector) length(x) == length(y) || throw(DimensionMismatch()) nzind = nonzeroinds(x) nzval = nonzeros(x) @@ -1481,7 +1481,7 @@ scale!(a::Complex, x::SparseVectorUnion) = (scale!(nonzeros(x), a); x) (/)(x::SparseVectorUnion, a::Number) = SparseVector(length(x), copy(nonzeroinds(x)), nonzeros(x) / a) # dot -function dot(x::StridedVector{Tx}, y::SparseVectorUnion{Ty}) where {Tx<:Number,Ty<:Number} +function dot(x::AbstractVector{Tx}, y::SparseVectorUnion{Ty}) where {Tx<:Number,Ty<:Number} n = length(x) length(y) == n || throw(DimensionMismatch()) nzind = nonzeroinds(y) @@ -1545,7 +1545,7 @@ end ### BLAS-2 / dense A * sparse x -> dense y # lowrankupdate (BLAS.ger! like) -function LinAlg.lowrankupdate!(A::StridedMatrix, x::StridedVector, y::SparseVectorUnion, α::Number = 1) +function LinAlg.lowrankupdate!(A::StridedMatrix, x::AbstractVector, y::SparseVectorUnion, α::Number = 1) nzi = nonzeroinds(y) nzv = nonzeros(y) @inbounds for (j,v) in zip(nzi,nzv) @@ -1567,10 +1567,10 @@ function (*)(A::StridedMatrix{Ta}, x::AbstractSparseVector{Tx}) where {Ta,Tx} mul!(y, A, x) end -mul!(y::StridedVector{Ty}, A::StridedMatrix, x::AbstractSparseVector{Tx}) where {Tx,Ty} = +mul!(y::AbstractVector{Ty}, A::StridedMatrix, x::AbstractSparseVector{Tx}) where {Tx,Ty} = mul!(one(Tx), A, x, zero(Ty), y) -function mul!(α::Number, A::StridedMatrix, x::AbstractSparseVector, β::Number, y::StridedVector) +function mul!(α::Number, A::StridedMatrix, x::AbstractSparseVector, β::Number, y::AbstractVector) m, n = size(A) length(x) == n && length(y) == m || throw(DimensionMismatch()) m == 0 && return y @@ -1605,10 +1605,10 @@ function *(transA::Transpose{<:Any,<:StridedMatrix{Ta}}, x::AbstractSparseVector mul!(y, Transpose(A), x) end -mul!(y::StridedVector{Ty}, transA::Transpose{<:Any,<:StridedMatrix}, x::AbstractSparseVector{Tx}) where {Tx,Ty} = +mul!(y::AbstractVector{Ty}, transA::Transpose{<:Any,<:StridedMatrix}, x::AbstractSparseVector{Tx}) where {Tx,Ty} = (A = transA.parent; mul!(one(Tx), Transpose(A), x, zero(Ty), y)) -function mul!(α::Number, transA::Transpose{<:Any,<:StridedMatrix}, x::AbstractSparseVector, β::Number, y::StridedVector) +function mul!(α::Number, transA::Transpose{<:Any,<:StridedMatrix}, x::AbstractSparseVector, β::Number, y::AbstractVector) A = transA.parent m, n = size(A) length(x) == m && length(y) == n || throw(DimensionMismatch()) @@ -1664,10 +1664,10 @@ end # * and mul! -mul!(y::StridedVector{Ty}, A::SparseMatrixCSC, x::AbstractSparseVector{Tx}) where {Tx,Ty} = +mul!(y::AbstractVector{Ty}, A::SparseMatrixCSC, x::AbstractSparseVector{Tx}) where {Tx,Ty} = mul!(one(Tx), A, x, zero(Ty), y) -function mul!(α::Number, A::SparseMatrixCSC, x::AbstractSparseVector, β::Number, y::StridedVector) +function mul!(α::Number, A::SparseMatrixCSC, x::AbstractSparseVector, β::Number, y::AbstractVector) m, n = size(A) length(x) == n && length(y) == m || throw(DimensionMismatch()) m == 0 && return y @@ -1697,21 +1697,21 @@ end # * and *(Tranpose(A), B) -mul!(y::StridedVector{Ty}, transA::Transpose{<:Any,<:SparseMatrixCSC}, x::AbstractSparseVector{Tx}) where {Tx,Ty} = +mul!(y::AbstractVector{Ty}, transA::Transpose{<:Any,<:SparseMatrixCSC}, x::AbstractSparseVector{Tx}) where {Tx,Ty} = (A = transA.parent; mul!(one(Tx), Transpose(A), x, zero(Ty), y)) -mul!(α::Number, transA::Transpose{<:Any,<:SparseMatrixCSC}, x::AbstractSparseVector, β::Number, y::StridedVector) = +mul!(α::Number, transA::Transpose{<:Any,<:SparseMatrixCSC}, x::AbstractSparseVector, β::Number, y::AbstractVector) = (A = transA.parent; _At_or_Ac_mul_B!(*, α, A, x, β, y)) -mul!(y::StridedVector{Ty}, adjA::Adjoint{<:Any,<:SparseMatrixCSC}, x::AbstractSparseVector{Tx}) where {Tx,Ty} = +mul!(y::AbstractVector{Ty}, adjA::Adjoint{<:Any,<:SparseMatrixCSC}, x::AbstractSparseVector{Tx}) where {Tx,Ty} = (A = adjA.parent; mul!(one(Tx), Adjoint(A), x, zero(Ty), y)) -mul!(α::Number, adjA::Adjoint{<:Any,<:SparseMatrixCSC}, x::AbstractSparseVector, β::Number, y::StridedVector) = +mul!(α::Number, adjA::Adjoint{<:Any,<:SparseMatrixCSC}, x::AbstractSparseVector, β::Number, y::AbstractVector) = (A = adjA.parent; _At_or_Ac_mul_B!(dot, α, A, x, β, y)) function _At_or_Ac_mul_B!(tfun::Function, α::Number, A::SparseMatrixCSC, x::AbstractSparseVector, - β::Number, y::StridedVector) + β::Number, y::AbstractVector) m, n = size(A) length(x) == m && length(y) == n || throw(DimensionMismatch()) n == 0 && return y diff --git a/base/subarray.jl b/base/subarray.jl index 6c8d79f72d726..d1a2b9bd94b1b 100644 --- a/base/subarray.jl +++ b/base/subarray.jl @@ -251,16 +251,16 @@ end IndexStyle(::Type{<:FastSubArray}) = IndexLinear() IndexStyle(::Type{<:SubArray}) = IndexCartesian() -# Strides are the distance between adjacent elements in a given dimension, -# so they are well-defined even for non-linear memory layouts +# Strides are the distance in memory between adjacent elements in a given dimension +# which we determine from the strides of the parent strides(V::SubArray) = substrides(V.parent, V.indices) -substrides(parent, I::Tuple) = substrides(1, parent, 1, I) -substrides(s, parent, dim, ::Tuple{}) = () -substrides(s, parent, dim, I::Tuple{ScalarIndex, Vararg{Any}}) = (substrides(s*size(parent, dim), parent, dim+1, tail(I))...,) -substrides(s, parent, dim, I::Tuple{Slice, Vararg{Any}}) = (s, substrides(s*size(parent, dim), parent, dim+1, tail(I))...) -substrides(s, parent, dim, I::Tuple{AbstractRange, Vararg{Any}}) = (s*step(I[1]), substrides(s*size(parent, dim), parent, dim+1, tail(I))...) -substrides(s, parent, dim, I::Tuple{Any, Vararg{Any}}) = throw(ArgumentError("strides is invalid for SubArrays with indices of type $(typeof(I[1]))")) +substrides(parent, I::Tuple) = substrides(parent, strides(parent), I) +substrides(parent, strds::Tuple{}, ::Tuple{}) = () +substrides(parent, strds::NTuple{N,Int}, I::Tuple{ScalarIndex, Vararg{Any}}) where N = (substrides(parent, tail(strds), tail(I))...,) +substrides(parent, strds::NTuple{N,Int}, I::Tuple{Slice, Vararg{Any}}) where N = (first(strds), substrides(parent, tail(strds), tail(I))...) +substrides(parent, strds::NTuple{N,Int}, I::Tuple{AbstractRange, Vararg{Any}}) where N = (first(strds)*step(I[1]), substrides(parent, tail(strds), tail(I))...) +substrides(parent, strds, I::Tuple{Any, Vararg{Any}}) = throw(ArgumentError("strides is invalid for SubArrays with indices of type $(typeof(I[1]))")) stride(V::SubArray, d::Integer) = d <= ndims(V) ? strides(V)[d] : strides(V)[end] * size(V)[end] diff --git a/doc/src/manual/interfaces.md b/doc/src/manual/interfaces.md index bf4d4ed09df07..c61104fbd2e4b 100644 --- a/doc/src/manual/interfaces.md +++ b/doc/src/manual/interfaces.md @@ -391,6 +391,37 @@ something other than 1), you should specialize `indices`. You should also specia so that the `dims` argument (ordinarily a `Dims` size-tuple) can accept `AbstractUnitRange` objects, perhaps range-types `Ind` of your own design. For more information, see [Arrays with custom indices](@ref). +## [Strided Arrays] + +| Methods to implement |   | Brief description | +|:----------------------------------------------- |:-------------------------------------- |:------------------------------------------------------------------------------------- | +| `strides(A)` |   | Return the distance in memory (in number of elements) between adjacent elements in each dimension as a tuple. If `A` is an `AbstractArray{T,0}`, this should return an empty tuple. | +| `Base.unsafe_convert(::Type{Ptr{T}}, A)` |   | Return the native address of an array. | +| **Optional methods** | **Default definition** | **Brief description** | +| `stride(A, i::Int)` |   `strides(A)[i]` | Return the distance in memory (in number of elements) between adjacent elements in dimension k. | + +A strided array is a subtype of `AbstractArray` whose entries are stored in memory with fixed strides. +Provided the element type of the array is compatible with BLAS, a strided array can utilize BLAS and LAPACK routines +for more efficient linear algebra routines. A typical example of a user-defined strided array is one +that wraps a standard `Array` with additional structure. + +Warning: do not implement these methods if the underlying storage is not actually strided, as it +may lead to incorrect results or segmentation faults. + +Here are some examples to demonstrate which type of arrays are strided and which are not: +```julia +1:5 # not strided (there is no storage associated with this array.) +Vector(1:5) # is strided with strides (1,) +A = [1 5; 2 6; 3 7; 4 8] # is strided with strides (1,4) +V = view(A, 1:2, :) # is strided with strides (1,4) +V = view(A, 1:2:3, 1:2) # is strided with strides (2,4) +V = view(A, [1,2,4], :) # is not strided, as the spacing between rows is not fixed. +``` + + + + + ## [Broadcasting](@id man-interfaces-broadcasting) | Methods to implement | Brief description | diff --git a/test/abstractarray.jl b/test/abstractarray.jl index a3149a9bb0223..2571bd617657c 100644 --- a/test/abstractarray.jl +++ b/test/abstractarray.jl @@ -431,13 +431,6 @@ function test_primitives(::Type{T}, shape, ::Type{TestAbstractArray}) where T # last(a) @test last(B) == B[length(B)] - # strides(a::AbstractArray) - @inferred strides(B) - strides_B = strides(B) - for (i, _stride) in enumerate(collect(strides_B)) - @test _stride == stride(B, i) - end - # isassigned(a::AbstractArray, i::Int...) j = rand(1:length(B)) @test isassigned(B, j) == true diff --git a/test/linalg/blas.jl b/test/linalg/blas.jl index 1c877e6464bcc..740bb944753d8 100644 --- a/test/linalg/blas.jl +++ b/test/linalg/blas.jl @@ -326,3 +326,112 @@ end @test triu(A[1,:] * A[1,:]') ≈ BLAS.her!('U', one(real(elty)), view(A, 1, :), zeros(elty, 5, 5)) @test tril(A[1,:] * A[1,:]') ≈ BLAS.her!('L', one(real(elty)), view(A, 1, :), zeros(elty, 5, 5)) end + +struct WrappedArray{T,N} <: AbstractArray{T,N} + A::Array{T,N} +end + +Base.size(A::WrappedArray) = size(A.A) +Base.getindex(A::WrappedArray, i::Int) = A.A[i] +Base.getindex(A::WrappedArray{T, N}, I::Vararg{Int, N}) where {T, N} = A.A[I...] +Base.setindex!(A::WrappedArray, v, i::Int) = setindex!(A.A, v, i) +Base.setindex!(A::WrappedArray{T, N}, v, I::Vararg{Int, N}) where {T, N} = setindex!(A.A, v, I...) +Base.unsafe_convert(::Type{Ptr{T}}, A::WrappedArray{T}) where T = Base.unsafe_convert(Ptr{T}, A.A) + +@test_deprecated strides(WrappedArray(rand(5))) +@test_deprecated stride(WrappedArray(rand(5)), 1) + +Base.stride(A::WrappedArray, i::Int) = stride(A.A, i) + +@testset "strided interface blas" begin + for elty in (Float32, Float64, ComplexF32, ComplexF64) + # Level 1 + x = WrappedArray(elty[1, 2, 3, 4]) + y = WrappedArray(elty[5, 6, 7, 8]) + BLAS.blascopy!(2, x, 1, y, 2) + @test y == WrappedArray(elty[1, 6, 2, 8]) + BLAS.scal!(2, elty(2), x, 1) + @test x == WrappedArray(elty[2, 4, 3, 4]) + @test BLAS.nrm2(1, x, 2) == elty(2) + @test BLAS.nrm2(x) == BLAS.nrm2(x.A) + BLAS.asum(x) == elty(13) + BLAS.axpy!(4, elty(2), x, 1, y, 1) + @test y == WrappedArray(elty[5, 14, 8, 16]) + BLAS.axpby!(elty(2), x, elty(3), y) + @test y == WrappedArray(elty[19, 50, 30, 56]) + @test BLAS.iamax(x) == 2 + # Level 2 + A = WrappedArray(elty[1 2; 3 4]) + x = WrappedArray(elty[1, 2]) + y = WrappedArray(elty[3, 4]) + @test BLAS.gemv!('N', elty(2), A, x, elty(1), y) isa WrappedArray{elty,1} + @test y == WrappedArray(elty[13, 26]) + @test BLAS.gbmv!('N', 2, 1, 0, elty(2), A, x, elty(1), y) isa WrappedArray{elty,1} + @test y == WrappedArray(elty[15, 40]) + @test BLAS.symv!('U', elty(2), A, x, elty(1), y) isa WrappedArray{elty,1} + @test y == WrappedArray(elty[25, 60]) + @test BLAS.trmv!('U', 'N', 'N', A, y) isa WrappedArray{elty,1} + @test y == WrappedArray(elty[145, 240]) + @test BLAS.trsv!('U', 'N', 'N', A, y) isa WrappedArray{elty,1} + @test y == WrappedArray(elty[25,60]) + @test BLAS.ger!(elty(2), x, y, A) isa WrappedArray{elty,2} + @test A == WrappedArray(elty[51 122; 103 244]) + @test BLAS.syr!('L', elty(2), x, A) isa WrappedArray{elty,2} + @test A == WrappedArray(elty[53 122; 107 252]) + # Level 3 + A = WrappedArray(elty[1 2; 3 4]) + B = WrappedArray(elty[5 6; 7 8]) + C = WrappedArray(elty[9 10; 11 12]) + BLAS.gemm!('N', 'N', elty(2), A, B, elty(1), C) isa WrappedArray{elty,2} + @test C == WrappedArray([47 54; 97 112]) + BLAS.symm!('L', 'U', elty(2), A, B, elty(1), C) isa WrappedArray{elty,2} + @test C == WrappedArray([85 98; 173 200]) + BLAS.syrk!('U', 'N', elty(2), A, elty(1), C) isa WrappedArray{elty,2} + @test C == WrappedArray([95 120; 173 250]) + BLAS.syr2k!('U', 'N', elty(2), A, B, elty(1), C) isa WrappedArray{elty,2} + @test C == WrappedArray([163 244; 173 462]) + BLAS.trmm!('L', 'U', 'N', 'N', elty(2), A, B) isa WrappedArray{elty,2} + @test B == WrappedArray([38 44; 56 64]) + BLAS.trsm!('L', 'U', 'N', 'N', elty(2), A, B) isa WrappedArray{elty,2} + @test B == WrappedArray([20 24; 28 32]) + end + for elty in (Float32, Float64) + # Level 1 + x = WrappedArray(elty[1, 2, 3, 4]) + y = WrappedArray(elty[5, 6, 7, 8]) + @test BLAS.dot(2, x, 1, y, 2) == elty(19) + # Level 2 + A = WrappedArray(elty[1 2; 3 4]) + x = WrappedArray(elty[1, 2]) + y = WrappedArray(elty[3, 4]) + BLAS.sbmv!('U', 1, elty(2), A, x, elty(1), y) isa WrappedArray{elty,1} + @test y == WrappedArray(elty[17,24]) + end + for elty in (ComplexF32, ComplexF64) + # Level 1 + x = WrappedArray(elty[1+im, 2+2im, 3+3im, 4+im]) + y = WrappedArray(elty[5-im, 6-2im, 7-3im, 8-im]) + @test BLAS.dotc(2, x, 1, y, 2) == elty(12-26im) + @test BLAS.dotu(2, x, 1, y, 2) == elty(26+12im) + # Level 2 + A = WrappedArray(elty[1+im 2+2im; 3+3im 4+4im]) + x = WrappedArray(elty[1+im, 2+2im]) + y = WrappedArray(elty[5-im, 6-2im]) + @test BLAS.hemv!('U', elty(2), A, x, elty(1), y) isa WrappedArray{elty,1} + @test y == WrappedArray(elty[7+17im, 30+14im]) + BLAS.hbmv!('U', 1, elty(2), A, x, elty(1), y) isa WrappedArray{elty,1} + @test y == WrappedArray(elty[13+39im, 54+30im]) + @test BLAS.her!('L', real(elty(2)), x, A) isa WrappedArray{elty,2} + @test A == WrappedArray(elty[5 2+2im; 11+3im 20]) + # Level 3 + A = WrappedArray(elty[1+im 2+2im; 3+3im 4+4im]) + B = WrappedArray(elty[1+im 2+2im; 3+3im 4+4im]) + C = WrappedArray(elty[1+im 2+2im; 3+3im 4+4im]) + @test BLAS.hemm!('L', 'U', elty(2), A, B, elty(1), C) isa WrappedArray{elty,2} + @test C == WrappedArray([3+27im 6+38im; 35+27im 52+36im]) + @test BLAS.herk!('U', 'N', real(elty(2)), A, real(elty(1)), C) isa WrappedArray{elty,2} + @test C == WrappedArray([23 50+38im; 35+27im 152]) + @test BLAS.her2k!('U', 'N', elty(2), A, B, real(elty(1)), C) isa WrappedArray{elty,2} + @test C == WrappedArray([63 138+38im; 35+27im 352]) + end +end diff --git a/test/linalg/dense.jl b/test/linalg/dense.jl index 55baf6ee2fae9..abf0a4f9f3cf7 100644 --- a/test/linalg/dense.jl +++ b/test/linalg/dense.jl @@ -817,9 +817,29 @@ end @test lyap(one(elty),a) == -a/2 end -@testset "stride1" begin +@testset "strides" begin a = rand(10) b = view(a,2:2:10) + A = rand(10,10) + B = view(A, 2:2:10, 2:2:10) + @test Base.LinAlg.stride1(a) == 1 @test Base.LinAlg.stride1(b) == 2 + + @test strides(a) == (1,) + @test strides(b) == (2,) + @test strides(A) == (1,10) + @test strides(B) == (2,20) + + @test_deprecated strides(1:5) + @test_deprecated stride(1:5,1) + + for M in (a, b, A, B) + @inferred strides(M) + strides_M = strides(M) + + for (i, _stride) in enumerate(collect(strides_M)) + @test _stride == stride(M, i) + end + end end diff --git a/test/linalg/lapack.jl b/test/linalg/lapack.jl index de22763d54d2d..2ec4d0bfa5762 100644 --- a/test/linalg/lapack.jl +++ b/test/linalg/lapack.jl @@ -4,6 +4,8 @@ using Test import Base.LinAlg.BlasInt + + @test_throws ArgumentError Base.LinAlg.LAPACK.chkuplo('Z') @test_throws ArgumentError Base.LinAlg.LAPACK.chkside('Z') @test_throws ArgumentError Base.LinAlg.LAPACK.chkdiag('Z') @@ -429,6 +431,11 @@ end @test_throws DimensionMismatch LAPACK.tgsen!(zeros(BlasInt,10),Z,Z,Z,zeros(elty,11,11)) @test_throws DimensionMismatch LAPACK.trsyl!('N','N',Z,Z,zeros(elty,11,11)) @test_throws DimensionMismatch LAPACK.tzrzf!(zeros(elty,10,5)) + + A = triu(rand(elty,4,4)) + V = view(A, 1:2, :) + M = Matrix(V) + @test LAPACK.tzrzf!(V) == LAPACK.tzrzf!(M) end end diff --git a/test/subarray.jl b/test/subarray.jl index ca84ae9347654..79a664044d246 100644 --- a/test/subarray.jl +++ b/test/subarray.jl @@ -347,7 +347,7 @@ x11289 = randn(5,5) # Tests where non-trailing dimensions are preserved A = copy(reshape(1:120, 3, 5, 8)) sA = view(A, 2:2, 1:5, :) -@test strides(sA) == (1, 3, 15) +@test @inferred(strides(sA)) == (1, 3, 15) @test parent(sA) == A @test parentindices(sA) == (2:2, 1:5, Base.Slice(1:8)) @test Base.parentdims(sA) == [1:3;] @@ -357,7 +357,7 @@ sA = view(A, 2:2, 1:5, :) sA[2:5:end] = -1 @test all(sA[2:5:end] .== -1) @test all(A[5:15:120] .== -1) -@test strides(sA) == (1,3,15) +@test @inferred(strides(sA)) == (1,3,15) @test stride(sA,3) == 15 @test stride(sA,4) == 120 test_bounds(sA) @@ -367,7 +367,7 @@ sA[1:3,1:5] = -2 @test all(A[:,:,5] .== -2) sA[:] = -3 @test all(A[:,:,5] .== -3) -@test strides(sA) == (1,3) +@test @inferred(strides(sA)) == (1,3) test_bounds(sA) sA = view(A, 1:3, 3:3, 2:5) @test Base.parentdims(sA) == [1:3;] @@ -378,7 +378,7 @@ sA = view(A, 1:3, 3:3, 2:5) test_bounds(sA) sA = view(A, 1:2:3, 1:3:5, 1:2:8) @test Base.parentdims(sA) == [1:3;] -@test strides(sA) == (2,9,30) +@test @inferred(strides(sA)) == (2,9,30) @test sA[:] == A[1:2:3, 1:3:5, 1:2:8][:] # issue #8807 @test view(view([1:5;], 1:5), 1:5) == [1:5;] @@ -409,7 +409,7 @@ sA = view(A, 2, :, 1:8) @test Base.parentdims(sA) == [2:3;] @test size(sA) == (5, 8) @test axes(sA) === (Base.OneTo(5), Base.OneTo(8)) -@test strides(sA) == (3,15) +@test @inferred(strides(sA)) == (3,15) @test sA[2, 1:8][:] == [5:15:120;] @test sA[:,1] == [2:3:14;] @test sA[2:5:end] == [5:15:110;] @@ -421,13 +421,13 @@ sA = view(A, 1:3, 1:5, 5) @test Base.parentdims(sA) == [1:2;] @test size(sA) == (3,5) @test axes(sA) === (Base.OneTo(3),Base.OneTo(5)) -@test strides(sA) == (1,3) +@test @inferred(strides(sA)) == (1,3) test_bounds(sA) sA = view(A, 1:2:3, 3, 1:2:8) @test Base.parentdims(sA) == [1,3] @test size(sA) == (2,4) @test axes(sA) === (Base.OneTo(2), Base.OneTo(4)) -@test strides(sA) == (2,30) +@test @inferred(strides(sA)) == (2,30) @test sA[:] == A[sA.indices...][:] test_bounds(sA) @@ -587,3 +587,10 @@ end @test sizeof(view(zeros(UInt8, 10), 1:4)) == 4 @test sizeof(view(zeros(UInt8, 10), 1:3)) == 3 @test sizeof(view(zeros(Float64, 10, 10), 1:3, 2:6)) == 120 + + +# PR #25321 +# checks that issue in type inference is resolved +A = rand(5,5,5,5) +V = view(A, 1:1 ,:, 1:3, :) +@test @inferred(strides(V)) == (1, 5, 25, 125)