diff --git a/NEWS.md b/NEWS.md index 61609b4539632..f44f4d3c03fd2 100644 --- a/NEWS.md +++ b/NEWS.md @@ -226,6 +226,10 @@ Language changes * `try` blocks without `catch` or `finally` are no longer allowed. An explicit empty `catch` block should be written instead ([#27554]). + * `AbstractArray` types that use unconventional (not 1-based) indexing can now support + `size`, `length`, and `@inbounds`. To optionally enforce conventional indices, + you can `@assert is_one_indexed(A)`. + Breaking changes ---------------- diff --git a/base/abstractarray.jl b/base/abstractarray.jl index 421a45724655e..f69addeb9cf46 100644 --- a/base/abstractarray.jl +++ b/base/abstractarray.jl @@ -83,6 +83,16 @@ function axes(A) map(OneTo, size(A)) end +""" + is_one_indexed(A) + is_one_indexed(A, B, ...) + +Return `true` if the indices of `A` start at 1 along each axis. +If multiple arrays are passed, all must be 1-indexed. +""" +is_one_indexed(A::AbstractArray) = all(x->first(x)==1, axes(A)) +is_one_indexed(A::AbstractArray, Bs::AbstractArray...) = all(is_one_indexed, (A, Bs...)) + # Performance optimization: get rid of a branch on `d` in `axes(A, d)` # for d=1. 1d arrays are heavily used, and the first dimension comes up # in other applications. diff --git a/base/array.jl b/base/array.jl index 77591d333d5bf..25eee1b8498c1 100644 --- a/base/array.jl +++ b/base/array.jl @@ -445,6 +445,7 @@ for (fname, felt) in ((:zeros, :zero), (:ones, :one)) end function _one(unit::T, x::AbstractMatrix) where T + @assert is_one_indexed(x) m,n = size(x) m==n || throw(DimensionMismatch("multiplicative identity defined only for square matrices")) # Matrix{T}(I, m, m) @@ -748,6 +749,7 @@ function setindex! end function setindex!(A::Array, X::AbstractArray, I::AbstractVector{Int}) @_propagate_inbounds_meta @boundscheck setindex_shape_check(X, length(I)) + @assert is_one_indexed(X) X′ = unalias(A, X) I′ = unalias(A, I) count = 1 @@ -876,6 +878,7 @@ append!(a::Vector, iter) = _append!(a, IteratorSize(iter), iter) push!(a::Vector, iter...) = append!(a, iter) function _append!(a, ::Union{HasLength,HasShape}, iter) + @assert is_one_indexed(a) n = length(a) resize!(a, n+length(iter)) @inbounds for (i,item) in zip(n+1:length(a), iter) @@ -923,6 +926,7 @@ prepend!(a::Vector, iter) = _prepend!(a, IteratorSize(iter), iter) pushfirst!(a::Vector, iter...) = prepend!(a, iter) function _prepend!(a, ::Union{HasLength,HasShape}, iter) + @assert is_one_indexed(a) n = length(iter) _growbeg!(a, n) i = 0 diff --git a/base/c.jl b/base/c.jl index a544a1ceaa120..1cd59b2ace603 100644 --- a/base/c.jl +++ b/base/c.jl @@ -267,6 +267,7 @@ transcode(T, src::String) = transcode(T, codeunits(src)) transcode(::Type{String}, src) = String(transcode(UInt8, src)) function transcode(::Type{UInt16}, src::AbstractVector{UInt8}) + @assert is_one_indexed(src) dst = UInt16[] i, n = 1, length(src) n > 0 || return dst @@ -317,6 +318,7 @@ function transcode(::Type{UInt16}, src::AbstractVector{UInt8}) end function transcode(::Type{UInt8}, src::AbstractVector{UInt16}) + @assert is_one_indexed(src) n = length(src) n == 0 && return UInt8[] diff --git a/base/combinatorics.jl b/base/combinatorics.jl index ed13bec2ac358..ac26f0f9f0ab6 100644 --- a/base/combinatorics.jl +++ b/base/combinatorics.jl @@ -64,6 +64,7 @@ isperm(p::Tuple{Int}) = p[1] == 1 isperm(p::Tuple{Int,Int}) = ((p[1] == 1) & (p[2] == 2)) | ((p[1] == 2) & (p[2] == 1)) function permute!!(a, p::AbstractVector{<:Integer}) + @assert is_one_indexed(a, p) count = 0 start = 0 while count < length(a) @@ -114,6 +115,7 @@ julia> A permute!(a, p::AbstractVector) = permute!!(a, copymutable(p)) function invpermute!!(a, p::AbstractVector{<:Integer}) + @assert is_one_indexed(a, p) count = 0 start = 0 while count < length(a) @@ -194,6 +196,7 @@ julia> B[invperm(v)] ``` """ function invperm(a::AbstractVector) + @assert is_one_indexed(a) b = zero(a) # similar vector of zeros n = length(a) @inbounds for (i, j) in enumerate(a) diff --git a/base/io.jl b/base/io.jl index 8b9ee4f0529ac..3a07fa2e0fb08 100644 --- a/base/io.jl +++ b/base/io.jl @@ -804,6 +804,7 @@ The size of `b` will be increased if needed (i.e. if `nb` is greater than `lengt and enough bytes could be read), but it will never be decreased. """ function readbytes!(s::IO, b::AbstractArray{UInt8}, nb=length(b)) + @assert is_one_indexed(b) olb = lb = length(b) nr = 0 while nr < nb && !eof(s) diff --git a/base/iobuffer.jl b/base/iobuffer.jl index c42093b9d2c59..ad50e3cc2ef28 100644 --- a/base/iobuffer.jl +++ b/base/iobuffer.jl @@ -16,6 +16,7 @@ mutable struct GenericIOBuffer{T<:AbstractVector{UInt8}} <: IO function GenericIOBuffer{T}(data::T, readable::Bool, writable::Bool, seekable::Bool, append::Bool, maxsize::Integer) where T<:AbstractVector{UInt8} + @assert is_one_indexed(data) new(data,readable,writable,seekable,append,length(data),maxsize,1,-1) end end @@ -167,6 +168,7 @@ function unsafe_read(from::GenericIOBuffer, p::Ptr{UInt8}, nb::UInt) end function read_sub(from::GenericIOBuffer, a::AbstractArray{T}, offs, nel) where T + @assert is_one_indexed(a) from.readable || throw(ArgumentError("read failed, IOBuffer is not readable")) if offs+nel-1 > length(a) || offs < 1 || nel < 0 throw(BoundsError()) @@ -401,6 +403,7 @@ function unsafe_write(to::GenericIOBuffer, p::Ptr{UInt8}, nb::UInt) end function write_sub(to::GenericIOBuffer, a::AbstractArray{UInt8}, offs, nel) + @assert is_one_indexed(a) if offs+nel-1 > length(a) || offs < 1 || nel < 0 throw(BoundsError()) end diff --git a/base/multidimensional.jl b/base/multidimensional.jl index 53c1f93d2b206..3dd8a79ba2515 100644 --- a/base/multidimensional.jl +++ b/base/multidimensional.jl @@ -647,7 +647,10 @@ _iterable(X::AbstractArray, I...) = X end end -diff(a::AbstractVector) = [ a[i+1] - a[i] for i=1:length(a)-1 ] +function diff(a::AbstractVector) + @assert is_one_indexed(a) + [ a[i+1] - a[i] for i=1:length(a)-1 ] +end """ diff(A::AbstractVector) @@ -1455,6 +1458,7 @@ function _extrema_dims(A::AbstractArray, dims) end @noinline function extrema!(B, A) + @assert is_one_indexed(B, A) sA = size(A) sB = size(B) for I in CartesianIndices(sB) diff --git a/base/number.jl b/base/number.jl index 8b819bd5739ed..1f75e51cc24c3 100644 --- a/base/number.jl +++ b/base/number.jl @@ -69,6 +69,7 @@ ndims(::Type{<:Number}) = 0 length(x::Number) = 1 firstindex(x::Number) = 1 lastindex(x::Number) = 1 +is_one_indexed(x::Number) = true IteratorSize(::Type{<:Number}) = HasShape{0}() keys(::Number) = OneTo(1) diff --git a/base/reshapedarray.jl b/base/reshapedarray.jl index f56b9268ddad0..57783b6cd8b3c 100644 --- a/base/reshapedarray.jl +++ b/base/reshapedarray.jl @@ -145,6 +145,7 @@ _reshape(parent::Array, dims::Dims) = reshape(parent, dims) # When reshaping Vector->Vector, don't wrap with a ReshapedArray function _reshape(v::AbstractVector, dims::Dims{1}) + @assert is_one_indexed(v) len = dims[1] len == length(v) || _throw_dmrs(_length(v), "length", len) v diff --git a/base/sort.jl b/base/sort.jl index 9f19b3b1e8b02..08bbc49cdb69b 100644 --- a/base/sort.jl +++ b/base/sort.jl @@ -225,6 +225,7 @@ function searchsorted(v::AbstractVector, x, ilo::Int, ihi::Int, o::Ordering) end function searchsortedlast(a::AbstractRange{<:Real}, x::Real, o::DirectOrdering) + is_one_indexed(a) || throw(ArgumentError("range must be indexed starting with 1")) if step(a) == 0 lt(o, x, first(a)) ? 0 : length(a) else @@ -234,6 +235,7 @@ function searchsortedlast(a::AbstractRange{<:Real}, x::Real, o::DirectOrdering) end function searchsortedfirst(a::AbstractRange{<:Real}, x::Real, o::DirectOrdering) + is_one_indexed(a) || throw(ArgumentError("range must be indexed starting with 1")) if step(a) == 0 lt(o, first(a), x) ? length(a) + 1 : 1 else @@ -243,6 +245,7 @@ function searchsortedfirst(a::AbstractRange{<:Real}, x::Real, o::DirectOrdering) end function searchsortedlast(a::AbstractRange{<:Integer}, x::Real, o::DirectOrdering) + is_one_indexed(a) || throw(ArgumentError("range must be indexed starting with 1")) if step(a) == 0 lt(o, x, first(a)) ? 0 : length(a) else @@ -251,6 +254,7 @@ function searchsortedlast(a::AbstractRange{<:Integer}, x::Real, o::DirectOrderin end function searchsortedfirst(a::AbstractRange{<:Integer}, x::Real, o::DirectOrdering) + is_one_indexed(a) || throw(ArgumentError("range must be indexed starting with 1")) if step(a) == 0 lt(o, first(a), x) ? length(a)+1 : 1 else @@ -259,6 +263,7 @@ function searchsortedfirst(a::AbstractRange{<:Integer}, x::Real, o::DirectOrderi end function searchsortedfirst(a::AbstractRange{<:Integer}, x::Unsigned, o::DirectOrdering) + is_one_indexed(a) || throw(ArgumentError("range must be indexed starting with 1")) if lt(o, first(a), x) if step(a) == 0 length(a) + 1 @@ -271,6 +276,7 @@ function searchsortedfirst(a::AbstractRange{<:Integer}, x::Unsigned, o::DirectOr end function searchsortedlast(a::AbstractRange{<:Integer}, x::Unsigned, o::DirectOrdering) + is_one_indexed(a) || throw(ArgumentError("range must be indexed starting with 1")) if lt(o, x, first(a)) 0 elseif step(a) == 0 diff --git a/base/tuple.jl b/base/tuple.jl index 0eade99391471..133ec7e6e57db 100644 --- a/base/tuple.jl +++ b/base/tuple.jl @@ -19,6 +19,7 @@ NTuple length(@nospecialize t::Tuple) = nfields(t) firstindex(@nospecialize t::Tuple) = 1 lastindex(@nospecialize t::Tuple) = length(t) +is_one_indexed(@nospecialize t::Tuple) = true size(@nospecialize(t::Tuple), d) = (d == 1) ? length(t) : throw(ArgumentError("invalid tuple dimension $d")) @eval getindex(t::Tuple, i::Int) = getfield(t, i, $(Expr(:boundscheck))) @eval getindex(t::Tuple, i::Real) = getfield(t, convert(Int, i), $(Expr(:boundscheck))) diff --git a/stdlib/Base64/src/decode.jl b/stdlib/Base64/src/decode.jl index 5e15de0bfbaca..89c4d03e14788 100644 --- a/stdlib/Base64/src/decode.jl +++ b/stdlib/Base64/src/decode.jl @@ -109,6 +109,7 @@ function Base.read(pipe::Base64DecodePipe, ::Type{UInt8}) end function Base.readbytes!(pipe::Base64DecodePipe, data::AbstractVector{UInt8}, nb::Integer=length(data)) + @assert is_one_indexed(data) filled::Int = 0 while filled < nb && !eof(pipe) if length(data) == filled diff --git a/stdlib/Distributed/src/pmap.jl b/stdlib/Distributed/src/pmap.jl index 9eb55dadd035a..d1839e8765e63 100644 --- a/stdlib/Distributed/src/pmap.jl +++ b/stdlib/Distributed/src/pmap.jl @@ -187,7 +187,7 @@ function wrap_batch(f, p, handle_errors) remotecall_fetch(f, p, batch) catch e if handle_errors - return Any[BatchProcessingError(batch[i], e) for i in 1:length(batch)] + return Any[BatchProcessingError(b, e) for b in batch] else rethrow(e) end diff --git a/stdlib/Future/src/Future.jl b/stdlib/Future/src/Future.jl index f87fca28bae6d..264c94c9abc36 100644 --- a/stdlib/Future/src/Future.jl +++ b/stdlib/Future/src/Future.jl @@ -23,8 +23,8 @@ copy!(dst::AbstractDict, src::AbstractDict) = merge!(empty!(dst), src) copy!(dst::AbstractVector, src::AbstractVector) = append!(empty!(dst), src) function copy!(dst::AbstractArray, src::AbstractArray) - size(dst) == size(src) || throw(ArgumentError( - "arrays must have the same size for copy! (consider using `copyto!`)")) + axes(dst) == axes(src) || throw(ArgumentError( + "arrays must have the same axes for copy! (consider using `copyto!`)")) copyto!(dst, src) end diff --git a/stdlib/LinearAlgebra/src/LinearAlgebra.jl b/stdlib/LinearAlgebra/src/LinearAlgebra.jl index d841f2fc8f2e5..902e2c4fe6a79 100644 --- a/stdlib/LinearAlgebra/src/LinearAlgebra.jl +++ b/stdlib/LinearAlgebra/src/LinearAlgebra.jl @@ -21,9 +21,6 @@ using Base: hvcat_fill, iszero, IndexLinear, _length, promote_op, promote_typeof @propagate_inbounds, @pure, reduce, typed_vcat using Base.Broadcast: Broadcasted -# We use `_length` because of non-1 indices; releases after julia 0.5 -# can go back to `length`. `_length(A)` is equivalent to `length(LinearIndices(A))`. - export # Modules LAPACK, diff --git a/stdlib/LinearAlgebra/src/adjtrans.jl b/stdlib/LinearAlgebra/src/adjtrans.jl index d4b08a07502a6..eef7aaf67b7a3 100644 --- a/stdlib/LinearAlgebra/src/adjtrans.jl +++ b/stdlib/LinearAlgebra/src/adjtrans.jl @@ -11,6 +11,7 @@ import Base: length, size, axes, IndexStyle, getindex, setindex!, parent, vec, c struct Adjoint{T,S} <: AbstractMatrix{T} parent::S function Adjoint{T,S}(A::S) where {T,S} + @assert is_one_indexed(A) checkeltype_adjoint(T, eltype(A)) new(A) end @@ -19,6 +20,7 @@ struct Transpose{T,S} <: AbstractMatrix{T} parent::S function Transpose{T,S}(A::S) where {T,S} checkeltype_transpose(T, eltype(A)) + @assert is_one_indexed(A) new(A) end end @@ -144,8 +146,8 @@ similar(A::AdjOrTransAbsVec) = wrapperop(A)(similar(A.parent)) similar(A::AdjOrTransAbsVec, ::Type{T}) where {T} = wrapperop(A)(similar(A.parent, Base.promote_op(wrapperop(A), T))) # for matrices, the semantics of the wrapped and unwrapped types are generally the same # and as you are allocating with similar anyway, you might as well get something unwrapped -similar(A::AdjOrTrans) = similar(A.parent, eltype(A), size(A)) -similar(A::AdjOrTrans, ::Type{T}) where {T} = similar(A.parent, T, size(A)) +similar(A::AdjOrTrans) = similar(A.parent, eltype(A), axes(A)) +similar(A::AdjOrTrans, ::Type{T}) where {T} = similar(A.parent, T, axes(A)) similar(A::AdjOrTrans, ::Type{T}, dims::Dims{N}) where {T,N} = similar(A.parent, T, dims) # sundry basic definitions @@ -197,6 +199,7 @@ broadcast(f, tvs::Union{Number,TransposeAbsVec}...) = transpose(broadcast((xs... *(u::AdjointAbsVec, v::AbstractVector) = dot(u.parent, v) *(u::TransposeAbsVec{T}, v::AbstractVector{T}) where {T<:Real} = dot(u.parent, v) function *(u::TransposeAbsVec, v::AbstractVector) + @assert is_one_indexed(u, v) @boundscheck length(u) == length(v) || throw(DimensionMismatch()) return sum(@inbounds(u[k]*v[k]) for k in 1:length(u)) end diff --git a/stdlib/LinearAlgebra/src/bidiag.jl b/stdlib/LinearAlgebra/src/bidiag.jl index ae48c605fada4..1cf5aed8c4558 100644 --- a/stdlib/LinearAlgebra/src/bidiag.jl +++ b/stdlib/LinearAlgebra/src/bidiag.jl @@ -6,6 +6,7 @@ struct Bidiagonal{T,V<:AbstractVector{T}} <: AbstractMatrix{T} ev::V # sub/super diagonal uplo::Char # upper bidiagonal ('U') or lower ('L') function Bidiagonal{T}(dv::V, ev::V, uplo::AbstractChar) where {T,V<:AbstractVector{T}} + @assert is_one_indexed(dv, ev) if length(ev) != length(dv)-1 throw(DimensionMismatch("length of diagonal vector is $(length(dv)), length of off-diagonal vector is $(length(ev))")) end @@ -352,6 +353,9 @@ mul!(C::AbstractMatrix, A::BiTri, B::Adjoint{<:Any,<:AbstractVecOrMat}) = A_mul_ mul!(C::AbstractVector, A::BiTri, B::Transpose{<:Any,<:AbstractVecOrMat}) = throw(MethodError(mul!, (C, A, B))) function check_A_mul_B!_sizes(C, A, B) + @assert is_one_indexed(C) + @assert is_one_indexed(A) + @assert is_one_indexed(B) nA, mA = size(A) nB, mB = size(B) nC, mC = size(C) @@ -432,6 +436,8 @@ function A_mul_B_td!(C::AbstractMatrix, A::BiTriSym, B::BiTriSym) end function A_mul_B_td!(C::AbstractVecOrMat, A::BiTriSym, B::AbstractVecOrMat) + @assert is_one_indexed(C) + @assert is_one_indexed(B) nA = size(A,1) nB = size(B,2) if !(size(C,1) == size(B,1) == nA) @@ -506,6 +512,7 @@ ldiv!(A::Union{Bidiagonal, AbstractTriangular}, b::AbstractVector) = naivesub!(A ldiv!(A::Transpose{<:Any,<:Bidiagonal}, b::AbstractVector) = ldiv!(copy(A), b) ldiv!(A::Adjoint{<:Any,<:Bidiagonal}, b::AbstractVector) = ldiv!(copy(A), b) function ldiv!(A::Union{Bidiagonal,AbstractTriangular}, B::AbstractMatrix) + @assert is_one_indexed(A, B) nA,mA = size(A) tmp = similar(B,size(B,1)) n = size(B, 1) @@ -520,6 +527,7 @@ function ldiv!(A::Union{Bidiagonal,AbstractTriangular}, B::AbstractMatrix) B end function ldiv!(adjA::Adjoint{<:Any,<:Union{Bidiagonal,AbstractTriangular}}, B::AbstractMatrix) + @assert is_one_indexed(adjA, B) A = adjA.parent nA,mA = size(A) tmp = similar(B,size(B,1)) @@ -535,6 +543,7 @@ function ldiv!(adjA::Adjoint{<:Any,<:Union{Bidiagonal,AbstractTriangular}}, B::A B end function ldiv!(transA::Transpose{<:Any,<:Union{Bidiagonal,AbstractTriangular}}, B::AbstractMatrix) + @assert is_one_indexed(transA, B) A = transA.parent nA,mA = size(A) tmp = similar(B,size(B,1)) @@ -551,6 +560,7 @@ function ldiv!(transA::Transpose{<:Any,<:Union{Bidiagonal,AbstractTriangular}}, end #Generic solver using naive substitution function naivesub!(A::Bidiagonal{T}, b::AbstractVector, x::AbstractVector = b) where T + @assert is_one_indexed(A, b, x) N = size(A, 2) if N != length(b) || N != length(x) throw(DimensionMismatch("second dimension of A, $N, does not match one of the lengths of x, $(length(x)), or b, $(length(b))")) diff --git a/stdlib/LinearAlgebra/src/blas.jl b/stdlib/LinearAlgebra/src/blas.jl index 90efaf1340633..7563153e128e9 100644 --- a/stdlib/LinearAlgebra/src/blas.jl +++ b/stdlib/LinearAlgebra/src/blas.jl @@ -321,6 +321,7 @@ for (fname, elty) in ((:cblas_zdotu_sub,:ComplexF64), end function dot(DX::Union{DenseArray{T},AbstractVector{T}}, DY::Union{DenseArray{T},AbstractVector{T}}) where T<:BlasReal + @assert is_one_indexed(DX, DY) n = length(DX) if n != length(DY) throw(DimensionMismatch("dot product arguments have lengths $(length(DX)) and $(length(DY))")) @@ -328,6 +329,7 @@ function dot(DX::Union{DenseArray{T},AbstractVector{T}}, DY::Union{DenseArray{T} GC.@preserve DX DY dot(n, pointer(DX), stride(DX, 1), pointer(DY), stride(DY, 1)) end function dotc(DX::Union{DenseArray{T},AbstractVector{T}}, DY::Union{DenseArray{T},AbstractVector{T}}) where T<:BlasComplex + @assert is_one_indexed(DX, DY) n = length(DX) if n != length(DY) throw(DimensionMismatch("dot product arguments have lengths $(length(DX)) and $(length(DY))")) @@ -335,6 +337,7 @@ function dotc(DX::Union{DenseArray{T},AbstractVector{T}}, DY::Union{DenseArray{T GC.@preserve DX DY dotc(n, pointer(DX), stride(DX, 1), pointer(DY), stride(DY, 1)) end function dotu(DX::Union{DenseArray{T},AbstractVector{T}}, DY::Union{DenseArray{T},AbstractVector{T}}) where T<:BlasComplex + @assert is_one_indexed(DX, DY) n = length(DX) if n != length(DY) throw(DimensionMismatch("dot product arguments have lengths $(length(DX)) and $(length(DY))")) @@ -515,6 +518,7 @@ for (fname, elty) in ((:daxpby_,:Float64), (:saxpby_,:Float32), end function axpby!(alpha::Number, x::Union{DenseArray{T},AbstractVector{T}}, beta::Number, y::Union{DenseArray{T},AbstractVector{T}}) where T<:BlasFloat + @assert is_one_indexed(x, y) if length(x) != length(y) throw(DimensionMismatch("x has length $(length(x)), but y has length $(length(y))")) end @@ -553,6 +557,7 @@ for (fname, elty) in ((:dgemv_,:Float64), #* .. Array Arguments .. # DOUBLE PRECISION A(LDA,*),X(*),Y(*) function gemv!(trans::AbstractChar, alpha::($elty), A::AbstractVecOrMat{$elty}, X::AbstractVector{$elty}, beta::($elty), Y::AbstractVector{$elty}) + @assert is_one_indexed(A, X, Y) 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))")) @@ -637,6 +642,7 @@ for (fname, elty) in ((:dgbmv_,:Float64), # * .. Array Arguments .. # DOUBLE PRECISION A(LDA,*),X(*),Y(*) function gbmv!(trans::AbstractChar, m::Integer, kl::Integer, ku::Integer, alpha::($elty), A::AbstractMatrix{$elty}, x::AbstractVector{$elty}, beta::($elty), y::AbstractVector{$elty}) + @assert is_one_indexed(A, x, y) chkstride1(A) ccall((@blasfunc($fname), libblas), Cvoid, (Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt}, Ref{BlasInt}, @@ -684,6 +690,7 @@ for (fname, elty, lib) in ((:dsymv_,:Float64,libblas), # .. Array Arguments .. # DOUBLE PRECISION A(LDA,*),X(*),Y(*) function symv!(uplo::AbstractChar, alpha::($elty), A::AbstractMatrix{$elty}, x::AbstractVector{$elty}, beta::($elty), y::AbstractVector{$elty}) + @assert is_one_indexed(A, x, y) m, n = size(A) if m != n throw(DimensionMismatch("matrix A is $m by $n but must be square")) @@ -735,6 +742,7 @@ for (fname, elty) in ((:zhemv_,:ComplexF64), (:chemv_,:ComplexF32)) @eval begin function hemv!(uplo::AbstractChar, α::$elty, A::AbstractMatrix{$elty}, x::AbstractVector{$elty}, β::$elty, y::AbstractVector{$elty}) + @assert is_one_indexed(A, x, y) m, n = size(A) if m != n throw(DimensionMismatch("matrix A is $m by $n but must be square")) @@ -779,6 +787,7 @@ for (fname, elty) in ((:dsbmv_,:Float64), # * .. Array Arguments .. # DOUBLE PRECISION A(LDA,*),X(*),Y(*) function sbmv!(uplo::AbstractChar, k::Integer, alpha::($elty), A::AbstractMatrix{$elty}, x::AbstractVector{$elty}, beta::($elty), y::AbstractVector{$elty}) + @assert is_one_indexed(A, x, y) chkstride1(A) ccall((@blasfunc($fname), libblas), Cvoid, (Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt}, Ref{$elty}, @@ -842,6 +851,7 @@ for (fname, elty) in ((:zhbmv_,:ComplexF64), # * .. Array Arguments .. # DOUBLE PRECISION A(LDA,*),X(*),Y(*) function hbmv!(uplo::AbstractChar, k::Integer, alpha::($elty), A::AbstractMatrix{$elty}, x::AbstractVector{$elty}, beta::($elty), y::AbstractVector{$elty}) + @assert is_one_indexed(A, x, y) chkstride1(A) ccall((@blasfunc($fname), libblas), Cvoid, (Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt}, Ref{$elty}, @@ -897,6 +907,7 @@ for (fname, elty) in ((:dtrmv_,:Float64), # * .. Array Arguments .. # DOUBLE PRECISION A(LDA,*),X(*) function trmv!(uplo::AbstractChar, trans::AbstractChar, diag::AbstractChar, A::AbstractMatrix{$elty}, x::AbstractVector{$elty}) + @assert is_one_indexed(A, B) n = checksquare(A) if n != length(x) throw(DimensionMismatch("A has size ($n,$n), x has length $(length(x))")) @@ -950,6 +961,7 @@ for (fname, elty) in ((:dtrsv_,:Float64), # .. Array Arguments .. # DOUBLE PRECISION A(LDA,*),X(*) function trsv!(uplo::AbstractChar, trans::AbstractChar, diag::AbstractChar, A::AbstractMatrix{$elty}, x::AbstractVector{$elty}) + @assert is_one_indexed(A, x) n = checksquare(A) if n != length(x) throw(DimensionMismatch("size of A is $n != length(x) = $(length(x))")) @@ -983,6 +995,7 @@ for (fname, elty) in ((:dger_,:Float64), (:cgerc_,:ComplexF32)) @eval begin function ger!(α::$elty, x::AbstractVector{$elty}, y::AbstractVector{$elty}, A::AbstractMatrix{$elty}) + @assert is_one_indexed(A, x, y) 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))")) @@ -1015,6 +1028,7 @@ for (fname, elty, lib) in ((:dsyr_,:Float64,libblas), (:csyr_,:ComplexF32,liblapack)) @eval begin function syr!(uplo::AbstractChar, α::$elty, x::AbstractVector{$elty}, A::AbstractMatrix{$elty}) + @assert is_one_indexed(A, x) n = checksquare(A) if length(x) != n throw(DimensionMismatch("A has size ($n,$n), x has length $(length(x))")) @@ -1044,6 +1058,7 @@ for (fname, elty, relty) in ((:zher_,:ComplexF64, :Float64), (:cher_,:ComplexF32, :Float32)) @eval begin function her!(uplo::AbstractChar, α::$relty, x::AbstractVector{$elty}, A::AbstractMatrix{$elty}) + @assert is_one_indexed(A, x) n = checksquare(A) if length(x) != n throw(DimensionMismatch("A has size ($n,$n), x has length $(length(x))")) @@ -1086,6 +1101,7 @@ for (gemm, elty) in # 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? + @assert is_one_indexed(A, B, C) m = size(A, transA == 'N' ? 1 : 2) ka = size(A, transA == 'N' ? 2 : 1) kb = size(B, transB == 'N' ? 1 : 2) @@ -1145,6 +1161,7 @@ for (mfname, elty) in ((:dsymm_,:Float64), # .. Array Arguments .. # DOUBLE PRECISION A(LDA,*),B(LDB,*),C(LDC,*) function symm!(side::AbstractChar, uplo::AbstractChar, alpha::($elty), A::AbstractMatrix{$elty}, B::AbstractMatrix{$elty}, beta::($elty), C::AbstractMatrix{$elty}) + @assert is_one_indexed(A, B, C) m, n = size(C) j = checksquare(A) if j != (side == 'L' ? m : n) @@ -1213,6 +1230,7 @@ for (mfname, elty) in ((:zhemm_,:ComplexF64), # .. Array Arguments .. # DOUBLE PRECISION A(LDA,*),B(LDB,*),C(LDC,*) function hemm!(side::AbstractChar, uplo::AbstractChar, alpha::($elty), A::AbstractMatrix{$elty}, B::AbstractMatrix{$elty}, beta::($elty), C::AbstractMatrix{$elty}) + @assert is_one_indexed(A, B, C) m, n = size(C) j = checksquare(A) if j != (side == 'L' ? m : n) @@ -1278,6 +1296,7 @@ for (fname, elty) in ((:dsyrk_,:Float64), function syrk!(uplo::AbstractChar, trans::AbstractChar, alpha::($elty), A::AbstractVecOrMat{$elty}, beta::($elty), C::AbstractMatrix{$elty}) + @assert is_one_indexed(A, C) 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 @@ -1334,6 +1353,7 @@ for (fname, elty, relty) in ((:zherk_, :ComplexF64, :Float64), # COMPLEX A(LDA,*),C(LDC,*) function herk!(uplo::AbstractChar, trans::AbstractChar, α::$relty, A::AbstractVecOrMat{$elty}, β::$relty, C::AbstractMatrix{$elty}) + @assert is_one_indexed(A, C) n = checksquare(C) nn = size(A, trans == 'N' ? 1 : 2) if nn != n @@ -1377,6 +1397,7 @@ for (fname, elty) in ((:dsyr2k_,:Float64), function syr2k!(uplo::AbstractChar, trans::AbstractChar, alpha::($elty), A::AbstractVecOrMat{$elty}, B::AbstractVecOrMat{$elty}, beta::($elty), C::AbstractMatrix{$elty}) + @assert is_one_indexed(A, B, C) 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 @@ -1417,6 +1438,7 @@ for (fname, elty1, elty2) in ((:zher2k_,:ComplexF64,:Float64), (:cher2k_,:Comple function her2k!(uplo::AbstractChar, trans::AbstractChar, alpha::($elty1), A::AbstractVecOrMat{$elty1}, B::AbstractVecOrMat{$elty1}, beta::($elty2), C::AbstractMatrix{$elty1}) + @assert is_one_indexed(A, B, C) 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 @@ -1504,6 +1526,7 @@ for (mmname, smname, elty) in # DOUBLE PRECISION A(LDA,*),B(LDB,*) function trmm!(side::AbstractChar, uplo::AbstractChar, transa::AbstractChar, diag::AbstractChar, alpha::Number, A::AbstractMatrix{$elty}, B::AbstractMatrix{$elty}) + @assert is_one_indexed(A, B, C) m, n = size(B) nA = checksquare(A) if nA != (side == 'L' ? m : n) @@ -1531,6 +1554,7 @@ for (mmname, smname, elty) in # DOUBLE PRECISION A(LDA,*),B(LDB,*) function trsm!(side::AbstractChar, uplo::AbstractChar, transa::AbstractChar, diag::AbstractChar, alpha::$elty, A::AbstractMatrix{$elty}, B::AbstractMatrix{$elty}) + @assert is_one_indexed(A, B) m, n = size(B) k = checksquare(A) if k != (side == 'L' ? m : n) diff --git a/stdlib/LinearAlgebra/src/bunchkaufman.jl b/stdlib/LinearAlgebra/src/bunchkaufman.jl index 19809dfa7c43d..a51cc862d7447 100644 --- a/stdlib/LinearAlgebra/src/bunchkaufman.jl +++ b/stdlib/LinearAlgebra/src/bunchkaufman.jl @@ -11,6 +11,11 @@ struct BunchKaufman{T,S<:AbstractMatrix} <: Factorization{T} symmetric::Bool rook::Bool info::BlasInt + + function BunchKaufman{T,S}(LD, ipiv, uplo, symmetric, rook, info) where {T,S<:AbstractMatrix} + @assert is_one_indexed(LD) + new(LD, ipiv, uplo, symmetric, rook, info) + end end BunchKaufman(A::AbstractMatrix{T}, ipiv::Vector{BlasInt}, uplo::AbstractChar, symmetric::Bool, rook::Bool, info::BlasInt) where {T} = @@ -121,6 +126,7 @@ issymmetric(B::BunchKaufman) = B.symmetric ishermitian(B::BunchKaufman) = !B.symmetric function _ipiv2perm_bk(v::AbstractVector{T}, maxi::Integer, uplo::AbstractChar) where T + @assert is_one_indexed(v) p = T[1:maxi;] uploL = uplo == 'L' i = uploL ? 1 : maxi diff --git a/stdlib/LinearAlgebra/src/cholesky.jl b/stdlib/LinearAlgebra/src/cholesky.jl index 8ab096276b3a1..f4d5ed9a1feec 100644 --- a/stdlib/LinearAlgebra/src/cholesky.jl +++ b/stdlib/LinearAlgebra/src/cholesky.jl @@ -32,6 +32,11 @@ struct Cholesky{T,S<:AbstractMatrix} <: Factorization{T} factors::S uplo::Char info::BlasInt + + function Cholesky{T,S}(factors, uplo, info) where {T,S<:AbstractMatrix} + @assert is_one_indexed(factors) + new(factors, uplo, info) + end end Cholesky(A::AbstractMatrix{T}, uplo::Symbol, info::BlasInt) where {T} = Cholesky{T,typeof(A)}(A, char_uplo(uplo), info) @@ -45,6 +50,11 @@ struct CholeskyPivoted{T,S<:AbstractMatrix} <: Factorization{T} rank::BlasInt tol::Real info::BlasInt + + function CholeskyPivoted{T,S}(factors, uplo, piv, rank, tol, info) where {T,S<:AbstractMatrix} + @assert is_one_indexed(factors) + new(factors, uplo, piv, rank, tol, info) + end end function CholeskyPivoted(A::AbstractMatrix{T}, uplo::AbstractChar, piv::Vector{BlasInt}, rank::BlasInt, tol::Real, info::BlasInt) where T @@ -75,6 +85,7 @@ end ## Non BLAS/LAPACK element types (generic) function _chol!(A::AbstractMatrix, ::Type{UpperTriangular}) + @assert is_one_indexed(A) n = checksquare(A) @inbounds begin for k = 1:n @@ -98,6 +109,7 @@ function _chol!(A::AbstractMatrix, ::Type{UpperTriangular}) return UpperTriangular(A), convert(BlasInt, 0) end function _chol!(A::AbstractMatrix, ::Type{LowerTriangular}) + @assert is_one_indexed(A) n = checksquare(A) @inbounds begin for k = 1:n diff --git a/stdlib/LinearAlgebra/src/dense.jl b/stdlib/LinearAlgebra/src/dense.jl index 0f422b70d44cb..df39a27d22351 100644 --- a/stdlib/LinearAlgebra/src/dense.jl +++ b/stdlib/LinearAlgebra/src/dense.jl @@ -169,6 +169,7 @@ julia> triu!(M, 1) ``` """ function triu!(M::AbstractMatrix, k::Integer) + @assert is_one_indexed(M) m, n = size(M) if !(-m + 1 <= k <= n + 1) throw(ArgumentError(string("the requested diagonal, $k, must be at least ", @@ -213,6 +214,7 @@ julia> tril!(M, 2) ``` """ function tril!(M::AbstractMatrix, k::Integer) + @assert is_one_indexed(M) m, n = size(M) if !(-m - 1 <= k <= n - 1) throw(ArgumentError(string("the requested diagonal, $k, must be at least ", @@ -236,6 +238,7 @@ tril(M::Matrix, k::Integer) = tril!(copy(M), k) Fill the band between diagonals `l` and `u` with the value `x`. """ function fillband!(A::AbstractMatrix{T}, x, l, u) where T + @assert is_one_indexed(A) m, n = size(A) xT = convert(T, x) for j in 1:n @@ -271,7 +274,10 @@ julia> diagind(A,-1) 2:4:6 ``` """ -diagind(A::AbstractMatrix, k::Integer=0) = diagind(size(A,1), size(A,2), k) +function diagind(A::AbstractMatrix, k::Integer=0) + @assert is_one_indexed(A) + diagind(size(A,1), size(A,2), k) +end """ diag(M, k::Integer=0) @@ -378,6 +384,7 @@ julia> kron(A, B) ``` """ function kron(a::AbstractMatrix{T}, b::AbstractMatrix{S}) where {T,S} + @assert is_one_indexed(a, b) R = Matrix{promote_op(*,T,S)}(undef, size(a,1)*size(b,1), size(a,2)*size(b,2)) m = 1 for j = 1:size(a,2), l = 1:size(b,2), i = 1:size(a,1) diff --git a/stdlib/LinearAlgebra/src/deprecated.jl b/stdlib/LinearAlgebra/src/deprecated.jl index 533a15d2f5fff..39ddae62bb2a5 100644 --- a/stdlib/LinearAlgebra/src/deprecated.jl +++ b/stdlib/LinearAlgebra/src/deprecated.jl @@ -186,6 +186,7 @@ _gradient(F::Vector, h::BitVector) = _gradient(F, Array(h)) _gradient(F::BitVector, h::Vector) = _gradient(Array(F), h) _gradient(F::BitVector, h::BitVector) = _gradient(Array(F), Array(h)) function _gradient(F::AbstractVector, h::Vector) + @assert is_one_indexed(F) n = length(F) T = typeof(oneunit(eltype(F))/oneunit(eltype(h))) g = similar(F, T) @@ -550,6 +551,7 @@ IndexStyle(::Type{<:RowVector}) = IndexLinear() # Generic behavior @inline function *(rowvec::RowVector, vec::AbstractVector) + @assert is_one_indexed(rowvec, vec) if length(rowvec) != length(vec) throw(DimensionMismatch("A has dimensions $(size(rowvec)) but B has dimensions $(size(vec))")) end diff --git a/stdlib/LinearAlgebra/src/diagonal.jl b/stdlib/LinearAlgebra/src/diagonal.jl index a6e528c6e827c..eac07b8964fc9 100644 --- a/stdlib/LinearAlgebra/src/diagonal.jl +++ b/stdlib/LinearAlgebra/src/diagonal.jl @@ -160,11 +160,13 @@ end lmul!(D, copyto!(similar(A, promote_op(*, eltype(A), eltype(D.diag)), size(A)), A)) function rmul!(A::AbstractMatrix, D::Diagonal) + @assert is_one_indexed(A) A .= A .* transpose(D.diag) return A end function lmul!(D::Diagonal, B::AbstractMatrix) + @assert is_one_indexed(B) B .= D.diag .* B return B end @@ -304,6 +306,7 @@ function ldiv!(D::Diagonal{T}, v::AbstractVector{T}) where {T} v end function ldiv!(D::Diagonal{T}, V::AbstractMatrix{T}) where {T} + @assert is_one_indexed(V) if size(V,1) != length(D.diag) throw(DimensionMismatch("diagonal matrix is $(length(D.diag)) by $(length(D.diag)) but right hand side has $(size(V,1)) rows")) end @@ -326,6 +329,7 @@ ldiv!(transD::Transpose{<:Any,<:Diagonal{T}}, B::AbstractVecOrMat{T}) where {T} (D = transD.parent; ldiv!(D, B)) function rdiv!(A::AbstractMatrix{T}, D::Diagonal{T}) where {T} + @assert is_one_indexed(A) dd = D.diag m, n = size(A) if (k = length(dd)) ≠ n diff --git a/stdlib/LinearAlgebra/src/factorization.jl b/stdlib/LinearAlgebra/src/factorization.jl index c41a55f65647c..f20653c3d2db2 100644 --- a/stdlib/LinearAlgebra/src/factorization.jl +++ b/stdlib/LinearAlgebra/src/factorization.jl @@ -73,18 +73,21 @@ end # With a real lhs and complex rhs with the same precision, we can reinterpret # the complex rhs as a real rhs with twice the number of columns function (\)(F::Factorization{T}, B::VecOrMat{Complex{T}}) where T<:BlasReal + @assert is_one_indexed(B) c2r = reshape(copy(transpose(reinterpret(T, reshape(B, (1, length(B)))))), size(B, 1), 2*size(B, 2)) x = ldiv!(F, c2r) return reshape(copy(reinterpret(Complex{T}, copy(transpose(reshape(x, div(length(x), 2), 2))))), _ret_size(F, B)) end function \(F::Factorization, B::AbstractVecOrMat) + @assert is_one_indexed(B) TFB = typeof(oneunit(eltype(B)) / oneunit(eltype(F))) BB = similar(B, TFB, size(B)) copyto!(BB, B) ldiv!(F, BB) end function \(adjF::Adjoint{<:Any,<:Factorization}, B::AbstractVecOrMat) + @assert is_one_indexed(B) F = adjF.parent TFB = typeof(oneunit(eltype(B)) / oneunit(eltype(F))) BB = similar(B, TFB, size(B)) @@ -94,6 +97,7 @@ end # support the same 3-arg idiom as in our other in-place A_*_B functions: function ldiv!(Y::AbstractVecOrMat, A::Factorization, B::AbstractVecOrMat) + @assert is_one_indexed(Y, B) m, n = size(A, 1), size(A, 2) if m > n ldiv!(A, B) diff --git a/stdlib/LinearAlgebra/src/generic.jl b/stdlib/LinearAlgebra/src/generic.jl index d0dc58607597a..61b994104bd7f 100644 --- a/stdlib/LinearAlgebra/src/generic.jl +++ b/stdlib/LinearAlgebra/src/generic.jl @@ -471,6 +471,7 @@ julia> norm(-2, Inf) # special cases of opnorm function opnorm1(A::AbstractMatrix{T}) where T + @assert is_one_indexed(A) m, n = size(A) Tnorm = typeof(float(real(zero(T)))) Tsum = promote_type(Float64, Tnorm) @@ -488,6 +489,7 @@ function opnorm1(A::AbstractMatrix{T}) where T end function opnorm2(A::AbstractMatrix{T}) where T + @assert is_one_indexed(A) m,n = size(A) if m == 1 || n == 1 return norm2(A) end Tnorm = typeof(float(real(zero(T)))) @@ -495,6 +497,7 @@ function opnorm2(A::AbstractMatrix{T}) where T end function opnormInf(A::AbstractMatrix{T}) where T + @assert is_one_indexed(A) m,n = size(A) Tnorm = typeof(float(real(zero(T)))) Tsum = promote_type(Float64, Tnorm) @@ -710,6 +713,7 @@ julia> dot([im; im], [1; 1]) ``` """ function dot(x::AbstractVector, y::AbstractVector) + @assert is_one_indexed(x, y) if length(LinearIndices(x)) != length(LinearIndices(y)) throw(DimensionMismatch("dot product arguments have unequal lengths $(length(LinearIndices(x))) and $(length(LinearIndices(y)))")) end @@ -872,6 +876,7 @@ true ``` """ function (\)(A::AbstractMatrix, B::AbstractVecOrMat) + @assert is_one_indexed(A, B) m, n = size(A) if m == n if istril(A) @@ -1031,6 +1036,7 @@ false ``` """ function istriu(A::AbstractMatrix, k::Integer = 0) + @assert is_one_indexed(A) m, n = size(A) for j in 1:min(n, m + k - 1) for i in max(1, j - k + 1):m @@ -1072,6 +1078,7 @@ false ``` """ function istril(A::AbstractMatrix, k::Integer = 0) + @assert is_one_indexed(A) m, n = size(A) for j in max(1, k + 2):n for i in 1:min(j - k - 1, m) @@ -1184,6 +1191,7 @@ end # Elementary reflection similar to LAPACK. The reflector is not Hermitian but # ensures that tridiagonalization of Hermitian matrices become real. See lawn72 @inline function reflector!(x::AbstractVector) + @assert is_one_indexed(x) n = length(x) @inbounds begin ξ1 = x[1] @@ -1207,6 +1215,7 @@ end # apply reflector from left @inline function reflectorApply!(x::AbstractVector, τ::Number, A::StridedMatrix) + @assert is_one_indexed(x) m, n = size(A) if length(x) != m throw(DimensionMismatch("reflector has length $(length(x)), which must match the first dimension of matrix A, $m")) diff --git a/stdlib/LinearAlgebra/src/givens.jl b/stdlib/LinearAlgebra/src/givens.jl index 49d45eb9e6f4e..412ee6fdc9e66 100644 --- a/stdlib/LinearAlgebra/src/givens.jl +++ b/stdlib/LinearAlgebra/src/givens.jl @@ -327,6 +327,7 @@ function getindex(G::Givens, i::Integer, j::Integer) end function lmul!(G::Givens, A::AbstractVecOrMat) + @assert is_one_indexed(A) m, n = size(A, 1), size(A, 2) if G.i2 > m throw(DimensionMismatch("column indices for rotation are outside the matrix")) @@ -339,6 +340,7 @@ function lmul!(G::Givens, A::AbstractVecOrMat) return A end function rmul!(A::AbstractMatrix, adjG::Adjoint{<:Any,<:Givens}) + @assert is_one_indexed(A) G = adjG.parent m, n = size(A, 1), size(A, 2) if G.i2 > n diff --git a/stdlib/LinearAlgebra/src/lapack.jl b/stdlib/LinearAlgebra/src/lapack.jl index 226c7aca5d222..63ea3469eff92 100644 --- a/stdlib/LinearAlgebra/src/lapack.jl +++ b/stdlib/LinearAlgebra/src/lapack.jl @@ -120,6 +120,7 @@ for (gbtrf, gbtrs, elty) in # INTEGER IPIV( * ) # DOUBLE PRECISION AB( LDAB, * ) function gbtrf!(kl::Integer, ku::Integer, m::Integer, AB::AbstractMatrix{$elty}) + @assert is_one_indexed(AB) chkstride1(AB) n = size(AB, 2) mnmn = min(m, n) @@ -143,6 +144,7 @@ for (gbtrf, gbtrs, elty) in function gbtrs!(trans::AbstractChar, kl::Integer, ku::Integer, m::Integer, AB::AbstractMatrix{$elty}, ipiv::AbstractVector{BlasInt}, B::AbstractVecOrMat{$elty}) + @assert is_one_indexed(AB, B) chkstride1(AB, B, ipiv) chktrans(trans) info = Ref{BlasInt}() @@ -221,6 +223,7 @@ for (gebal, gebak, elty, relty) in function gebak!(job::AbstractChar, side::AbstractChar, ilo::BlasInt, ihi::BlasInt, scale::AbstractVector{$relty}, V::AbstractMatrix{$elty}) + @assert is_one_indexed(scale, V) chkstride1(scale, V) chkside(side) chkfinite(V) # balancing routines don't support NaNs and Infs @@ -279,6 +282,7 @@ for (gebrd, gelqf, geqlf, geqrf, geqp3, geqrt, geqrt3, gerqf, getrf, elty, relty # DOUBLE PRECISION A( LDA, * ), D( * ), E( * ), TAUP( * ), # TAUQ( * ), WORK( * ) function gebrd!(A::AbstractMatrix{$elty}) + @assert is_one_indexed(A) chkstride1(A) m, n = size(A) k = min(m, n) @@ -312,6 +316,7 @@ for (gebrd, gelqf, geqlf, geqrf, geqp3, geqrt, geqrt3, gerqf, getrf, elty, relty # * .. Array Arguments .. # DOUBLE PRECISION A( LDA, * ), TAU( * ), WORK( * ) function gelqf!(A::AbstractMatrix{$elty}, tau::AbstractVector{$elty}) + @assert is_one_indexed(A, tau) chkstride1(A,tau) m = BlasInt(size(A, 1)) n = BlasInt(size(A, 2)) @@ -342,6 +347,7 @@ for (gebrd, gelqf, geqlf, geqrf, geqp3, geqrt, geqrt3, gerqf, getrf, elty, relty # * .. Array Arguments .. # DOUBLE PRECISION A( LDA, * ), TAU( * ), WORK( * ) function geqlf!(A::AbstractMatrix{$elty}, tau::AbstractVector{$elty}) + @assert is_one_indexed(A, tau) chkstride1(A,tau) m = BlasInt(size(A, 1)) n = BlasInt(size(A, 2)) @@ -373,6 +379,7 @@ for (gebrd, gelqf, geqlf, geqrf, geqp3, geqrt, geqrt3, gerqf, getrf, elty, relty # INTEGER JPVT( * ) # DOUBLE PRECISION A( LDA, * ), TAU( * ), WORK( * ) function geqp3!(A::AbstractMatrix{$elty}, jpvt::AbstractVector{BlasInt}, tau::AbstractVector{$elty}) + @assert is_one_indexed(A, jpvt, tau) chkstride1(A,jpvt,tau) m,n = size(A) if length(tau) != min(m,n) @@ -420,6 +427,7 @@ for (gebrd, gelqf, geqlf, geqrf, geqp3, geqrt, geqrt3, gerqf, getrf, elty, relty end function geqrt!(A::AbstractMatrix{$elty}, T::AbstractMatrix{$elty}) + @assert is_one_indexed(A, T) chkstride1(A) m, n = size(A) minmn = min(m, n) @@ -444,6 +452,7 @@ for (gebrd, gelqf, geqlf, geqrf, geqp3, geqrt, geqrt3, gerqf, getrf, elty, relty end function geqrt3!(A::AbstractMatrix{$elty}, T::AbstractMatrix{$elty}) + @assert is_one_indexed(A, T) chkstride1(A) chkstride1(T) m, n = size(A) @@ -473,6 +482,7 @@ for (gebrd, gelqf, geqlf, geqrf, geqp3, geqrt, geqrt3, gerqf, getrf, elty, relty # * .. Array Arguments .. # DOUBLE PRECISION A( LDA, * ), TAU( * ), WORK( * ) function geqrf!(A::AbstractMatrix{$elty}, tau::AbstractVector{$elty}) + @assert is_one_indexed(A, tau) chkstride1(A,tau) m, n = size(A) if length(tau) != min(m,n) @@ -501,6 +511,7 @@ for (gebrd, gelqf, geqlf, geqrf, geqp3, geqrt, geqrt3, gerqf, getrf, elty, relty # * .. Array Arguments .. # DOUBLE PRECISION A( LDA, * ), TAU( * ), WORK( * ) function gerqf!(A::AbstractMatrix{$elty},tau::AbstractVector{$elty}) + @assert is_one_indexed(A, tau) chkstride1(A,tau) m, n = size(A) if length(tau) != min(m,n) @@ -530,6 +541,7 @@ for (gebrd, gelqf, geqlf, geqrf, geqp3, geqrt, geqrt3, gerqf, getrf, elty, relty # INTEGER IPIV( * ) # DOUBLE PRECISION A( LDA, * ) function getrf!(A::AbstractMatrix{$elty}) + @assert is_one_indexed(A) chkstride1(A) m, n = size(A) lda = max(1,stride(A, 2)) @@ -756,6 +768,7 @@ for (tzrzf, ormrz, elty) in # .. Array Arguments .. # COMPLEX*16 A( LDA, * ), TAU( * ), WORK( * ) function tzrzf!(A::AbstractMatrix{$elty}) + @assert is_one_indexed(A) chkstride1(A) m, n = size(A) if n < m @@ -792,6 +805,7 @@ for (tzrzf, ormrz, elty) in # COMPLEX*16 A( LDA, * ), C( LDC, * ), TAU( * ), WORK( * ) function ormrz!(side::AbstractChar, trans::AbstractChar, A::AbstractMatrix{$elty}, tau::AbstractVector{$elty}, C::AbstractMatrix{$elty}) + @assert is_one_indexed(A, tau, C) chktrans(trans) chkside(side) chkstride1(A, tau, C) @@ -857,6 +871,7 @@ for (gels, gesv, getrs, getri, elty) in # CHARACTER TRANS # INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS function gels!(trans::AbstractChar, A::AbstractMatrix{$elty}, B::AbstractVecOrMat{$elty}) + @assert is_one_indexed(A, B) chktrans(trans) chkstride1(A, B) btrn = trans == 'T' @@ -901,6 +916,7 @@ for (gels, gesv, getrs, getri, elty) in # INTEGER IPIV( * ) # DOUBLE PRECISION A( LDA, * ), B( LDB, * ) function gesv!(A::AbstractMatrix{$elty}, B::AbstractVecOrMat{$elty}) + @assert is_one_indexed(A, B) chkstride1(A, B) n = checksquare(A) if size(B,1) != n @@ -924,6 +940,7 @@ for (gels, gesv, getrs, getri, elty) in # INTEGER IPIV( * ) # DOUBLE PRECISION A( LDA, * ), B( LDB, * ) function getrs!(trans::AbstractChar, A::AbstractMatrix{$elty}, ipiv::AbstractVector{BlasInt}, B::AbstractVecOrMat{$elty}) + @assert is_one_indexed(A, ipiv, B) chktrans(trans) chkstride1(A, B, ipiv) n = checksquare(A) @@ -947,6 +964,7 @@ for (gels, gesv, getrs, getri, elty) in # INTEGER IPIV( * ) # DOUBLE PRECISION A( LDA, * ), WORK( * ) function getri!(A::AbstractMatrix{$elty}, ipiv::AbstractVector{BlasInt}) + @assert is_one_indexed(A, ipiv) chkstride1(A, ipiv) n = checksquare(A) if n != length(ipiv) @@ -1037,6 +1055,7 @@ for (gesvx, elty) in function gesvx!(fact::AbstractChar, trans::AbstractChar, A::AbstractMatrix{$elty}, AF::AbstractMatrix{$elty}, ipiv::AbstractVector{BlasInt}, equed::AbstractChar, R::AbstractVector{$elty}, C::AbstractVector{$elty}, B::AbstractVecOrMat{$elty}) + @assert is_one_indexed(A, AF, ipiv, R, C, B) chktrans(trans) chkstride1(ipiv, R, C, B) n = checksquare(A) @@ -1106,6 +1125,7 @@ for (gesvx, elty, relty) in function gesvx!(fact::AbstractChar, trans::AbstractChar, A::AbstractMatrix{$elty}, AF::AbstractMatrix{$elty}, ipiv::AbstractVector{BlasInt}, equed::AbstractChar, R::AbstractVector{$relty}, C::AbstractVector{$relty}, B::AbstractVecOrMat{$elty}) + @assert is_one_indexed(A, AF, ipiv, R, C, B) chktrans(trans) chkstride1(A, AF, ipiv, R, C, B) n = checksquare(A) @@ -1203,6 +1223,7 @@ for (gelsd, gelsy, elty) in # INTEGER IWORK( * ) # DOUBLE PRECISION A( LDA, * ), B( LDB, * ), S( * ), WORK( * ) function gelsd!(A::AbstractMatrix{$elty}, B::AbstractVecOrMat{$elty}, rcond::Real=-one($elty)) + @assert is_one_indexed(A, B) chkstride1(A, B) m, n = size(A) if size(B, 1) != m @@ -1245,6 +1266,7 @@ for (gelsd, gelsy, elty) in # INTEGER JPVT( * ) # DOUBLE PRECISION A( LDA, * ), B( LDB, * ), WORK( * ) function gelsy!(A::AbstractMatrix{$elty}, B::AbstractVecOrMat{$elty}, rcond::Real=eps($elty)) + @assert is_one_indexed(A, B) chkstride1(A) m = size(A, 1) n = size(A, 2) @@ -1296,6 +1318,7 @@ for (gelsd, gelsy, elty, relty) in # DOUBLE PRECISION RWORK( * ), S( * ) # COMPLEX*16 A( LDA, * ), B( LDB, * ), WORK( * ) function gelsd!(A::AbstractMatrix{$elty}, B::AbstractVecOrMat{$elty}, rcond::Real=-one($relty)) + @assert is_one_indexed(A, B) chkstride1(A, B) m, n = size(A) if size(B, 1) != m @@ -1341,6 +1364,7 @@ for (gelsd, gelsy, elty, relty) in # DOUBLE PRECISION RWORK( * ) # COMPLEX*16 A( LDA, * ), B( LDB, * ), WORK( * ) function gelsy!(A::AbstractMatrix{$elty}, B::AbstractVecOrMat{$elty}, rcond::Real=eps($relty)) + @assert is_one_indexed(A, B) chkstride1(A, B) m, n = size(A) nrhs = size(B, 2) @@ -1414,6 +1438,7 @@ for (gglse, elty) in ((:dgglse_, :Float64), # $ WORK( * ), X( * ) function gglse!(A::AbstractMatrix{$elty}, c::AbstractVector{$elty}, B::AbstractMatrix{$elty}, d::AbstractVector{$elty}) + @assert is_one_indexed(A, c, B, d) chkstride1(A, c, B, d) m, n = size(A) p = size(B, 1) @@ -1530,6 +1555,7 @@ for (geev, gesvd, gesdd, ggsvd, elty, relty) in # DOUBLE PRECISION A( LDA, * ), S( * ), U( LDU, * ), # VT( LDVT, * ), WORK( * ) function gesdd!(job::AbstractChar, A::AbstractMatrix{$elty}) + @assert is_one_indexed(A) chkstride1(A) m, n = size(A) minmn = min(m, n) @@ -1610,6 +1636,7 @@ for (geev, gesvd, gesdd, ggsvd, elty, relty) in # DOUBLE PRECISION A( LDA, * ), S( * ), U( LDU, * ), # $ VT( LDVT, * ), WORK( * ) function gesvd!(jobu::AbstractChar, jobvt::AbstractChar, A::AbstractMatrix{$elty}) + @assert is_one_indexed(A) chkstride1(A) m, n = size(A) minmn = min(m, n) @@ -1678,6 +1705,7 @@ for (geev, gesvd, gesdd, ggsvd, elty, relty) in # COMPLEX*16 A( LDA, * ), B( LDB, * ), Q( LDQ, * ), # $ U( LDU, * ), V( LDV, * ), WORK( * ) function ggsvd!(jobu::AbstractChar, jobv::AbstractChar, jobq::AbstractChar, A::AbstractMatrix{$elty}, B::AbstractMatrix{$elty}) + @assert is_one_indexed(A, B) chkstride1(A, B) m, n = size(A) if size(B, 2) != n @@ -1800,6 +1828,7 @@ for (f, elty) in ((:dggsvd3_, :Float64), (:sggsvd3_, :Float32)) @eval begin function ggsvd3!(jobu::AbstractChar, jobv::AbstractChar, jobq::AbstractChar, A::AbstractMatrix{$elty}, B::AbstractMatrix{$elty}) + @assert is_one_indexed(A, B) chkstride1(A, B) m, n = size(A) if size(B, 2) != n @@ -1856,6 +1885,7 @@ for (f, elty, relty) in ((:zggsvd3_, :ComplexF64, :Float64), (:cggsvd3_, :ComplexF32, :Float32)) @eval begin function ggsvd3!(jobu::AbstractChar, jobv::AbstractChar, jobq::AbstractChar, A::AbstractMatrix{$elty}, B::AbstractMatrix{$elty}) + @assert is_one_indexed(A, B) chkstride1(A, B) m, n = size(A) if size(B, 2) != n @@ -2021,6 +2051,7 @@ for (geevx, ggev, elty) in # $ B( LDB, * ), BETA( * ), VL( LDVL, * ), # $ VR( LDVR, * ), WORK( * ) function ggev!(jobvl::AbstractChar, jobvr::AbstractChar, A::AbstractMatrix{$elty}, B::AbstractMatrix{$elty}) + @assert is_one_indexed(A, B) chkstride1(A,B) n, m = checksquare(A,B) if n != m @@ -2167,6 +2198,7 @@ for (geevx, ggev, elty, relty) in # $ BETA( * ), VL( LDVL, * ), VR( LDVR, * ), # $ WORK( * ) function ggev!(jobvl::AbstractChar, jobvr::AbstractChar, A::AbstractMatrix{$elty}, B::AbstractMatrix{$elty}) + @assert is_one_indexed(A, B) chkstride1(A, B) n, m = checksquare(A, B) if n != m @@ -2265,6 +2297,7 @@ for (laic1, elty) in # DOUBLE PRECISION W( J ), X( J ) function laic1!(job::Integer, x::AbstractVector{$elty}, sest::$elty, w::AbstractVector{$elty}, gamma::$elty) + @assert is_one_indexed(x, w) 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))")) @@ -2298,6 +2331,7 @@ for (laic1, elty, relty) in # COMPLEX*16 W( J ), X( J ) function laic1!(job::Integer, x::AbstractVector{$elty}, sest::$relty, w::AbstractVector{$elty}, gamma::$elty) + @assert is_one_indexed(x, w) 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))")) @@ -2331,6 +2365,7 @@ for (gtsv, gttrf, gttrs, elty) in # DOUBLE PRECISION B( LDB, * ), D( * ), DL( * ), DU( * ) function gtsv!(dl::AbstractVector{$elty}, d::AbstractVector{$elty}, du::AbstractVector{$elty}, B::AbstractVecOrMat{$elty}) + @assert is_one_indexed(dl, d, du, B) chkstride1(B, dl, d, du) n = length(d) if !(n >= length(dl) >= n - 1) @@ -2361,6 +2396,7 @@ for (gtsv, gttrf, gttrs, elty) in # INTEGER IPIV( * ) # DOUBLE PRECISION D( * ), DL( * ), DU( * ), DU2( * ) function gttrf!(dl::AbstractVector{$elty}, d::AbstractVector{$elty}, du::AbstractVector{$elty}) + @assert is_one_indexed(dl, d, du) chkstride1(dl,d,du) n = length(d) if length(dl) != n - 1 @@ -2390,6 +2426,7 @@ for (gtsv, gttrf, gttrs, elty) in function gttrs!(trans::AbstractChar, dl::AbstractVector{$elty}, d::AbstractVector{$elty}, du::AbstractVector{$elty}, du2::AbstractVector{$elty}, ipiv::AbstractVector{BlasInt}, B::AbstractVecOrMat{$elty}) + @assert is_one_indexed(dl, d, du, du2, ipiv, B) chktrans(trans) chkstride1(B, ipiv, dl, d, du, du2) n = length(d) @@ -2459,6 +2496,7 @@ for (orglq, orgqr, orgql, orgrq, ormlq, ormqr, ormql, ormrq, gemqrt, elty) in # * .. Array Arguments .. # DOUBLE PRECISION A( LDA, * ), TAU( * ), WORK( * ) function orglq!(A::AbstractMatrix{$elty}, tau::AbstractVector{$elty}, k::Integer = length(tau)) + @assert is_one_indexed(A, tau) chkstride1(A,tau) n = size(A, 2) m = min(n, size(A, 1)) @@ -2492,6 +2530,7 @@ for (orglq, orgqr, orgql, orgrq, ormlq, ormqr, ormql, ormrq, gemqrt, elty) in # * .. Array Arguments .. # DOUBLE PRECISION A( LDA, * ), TAU( * ), WORK( * ) function orgqr!(A::AbstractMatrix{$elty}, tau::AbstractVector{$elty}, k::Integer = length(tau)) + @assert is_one_indexed(A, tau) chkstride1(A,tau) m = size(A, 1) n = min(m, size(A, 2)) @@ -2527,6 +2566,7 @@ for (orglq, orgqr, orgql, orgrq, ormlq, ormqr, ormql, ormrq, gemqrt, elty) in # * .. Array Arguments .. # DOUBLE PRECISION A( LDA, * ), TAU( * ), WORK( * ) function orgql!(A::AbstractMatrix{$elty}, tau::AbstractVector{$elty}, k::Integer = length(tau)) + @assert is_one_indexed(A, tau) chkstride1(A,tau) m = size(A, 1) n = min(m, size(A, 2)) @@ -2562,6 +2602,7 @@ for (orglq, orgqr, orgql, orgrq, ormlq, ormqr, ormql, ormrq, gemqrt, elty) in # * .. Array Arguments .. # DOUBLE PRECISION A( LDA, * ), TAU( * ), WORK( * ) function orgrq!(A::AbstractMatrix{$elty}, tau::AbstractVector{$elty}, k::Integer = length(tau)) + @assert is_one_indexed(A, tau) chkstride1(A,tau) m, n = size(A) if n < m @@ -2598,6 +2639,7 @@ for (orglq, orgqr, orgql, orgrq, ormlq, ormqr, ormql, ormrq, gemqrt, elty) in # DOUBLE PRECISION A( LDA, * ), C( LDC, * ), TAU( * ), WORK( * ) function ormlq!(side::AbstractChar, trans::AbstractChar, A::AbstractMatrix{$elty}, tau::AbstractVector{$elty}, C::AbstractVecOrMat{$elty}) + @assert is_one_indexed(A, tau, C) chktrans(trans) chkside(side) chkstride1(A, C, tau) @@ -2644,6 +2686,7 @@ for (orglq, orgqr, orgql, orgrq, ormlq, ormqr, ormql, ormrq, gemqrt, elty) in # DOUBLE PRECISION A( LDA, * ), C( LDC, * ), TAU( * ), WORK( * ) function ormqr!(side::AbstractChar, trans::AbstractChar, A::AbstractMatrix{$elty}, tau::AbstractVector{$elty}, C::AbstractVecOrMat{$elty}) + @assert is_one_indexed(A, tau, C) chktrans(trans) chkside(side) chkstride1(A, C, tau) @@ -2693,6 +2736,7 @@ for (orglq, orgqr, orgql, orgrq, ormlq, ormqr, ormql, ormrq, gemqrt, elty) in # DOUBLE PRECISION A( LDA, * ), C( LDC, * ), TAU( * ), WORK( * ) function ormql!(side::AbstractChar, trans::AbstractChar, A::AbstractMatrix{$elty}, tau::AbstractVector{$elty}, C::AbstractVecOrMat{$elty}) + @assert is_one_indexed(A, tau, C) chktrans(trans) chkside(side) chkstride1(A, C, tau) @@ -2742,6 +2786,7 @@ for (orglq, orgqr, orgql, orgrq, ormlq, ormqr, ormql, ormrq, gemqrt, elty) in # DOUBLE PRECISION A( LDA, * ), C( LDC, * ), TAU( * ), WORK( * ) function ormrq!(side::AbstractChar, trans::AbstractChar, A::AbstractMatrix{$elty}, tau::AbstractVector{$elty}, C::AbstractVecOrMat{$elty}) + @assert is_one_indexed(A, tau, C) chktrans(trans) chkside(side) chkstride1(A, C, tau) @@ -2780,6 +2825,7 @@ for (orglq, orgqr, orgql, orgrq, ormlq, ormqr, ormql, ormrq, gemqrt, elty) in end function gemqrt!(side::AbstractChar, trans::AbstractChar, V::AbstractMatrix{$elty}, T::AbstractMatrix{$elty}, C::AbstractVecOrMat{$elty}) + @assert is_one_indexed(V, T, C) chktrans(trans) chkside(side) chkstride1(V, T, C) @@ -2930,6 +2976,7 @@ for (posv, potrf, potri, potrs, pstrf, elty, rtyp) in # .. Array Arguments .. # DOUBLE PRECISION A( LDA, * ), B( LDB, * ) function posv!(uplo::AbstractChar, A::AbstractMatrix{$elty}, B::AbstractVecOrMat{$elty}) + @assert is_one_indexed(A, B) chkstride1(A, B) n = checksquare(A) chkuplo(uplo) @@ -2953,6 +3000,7 @@ for (posv, potrf, potri, potrs, pstrf, elty, rtyp) in # * .. Array Arguments .. # DOUBLE PRECISION A( LDA, * ) function potrf!(uplo::AbstractChar, A::AbstractMatrix{$elty}) + @assert is_one_indexed(A) chkstride1(A) checksquare(A) chkuplo(uplo) @@ -2978,6 +3026,7 @@ for (posv, potrf, potri, potrs, pstrf, elty, rtyp) in # .. Array Arguments .. # DOUBLE PRECISION A( LDA, * ) function potri!(uplo::AbstractChar, A::AbstractMatrix{$elty}) + @assert is_one_indexed(A) chkstride1(A) chkuplo(uplo) info = Ref{BlasInt}() @@ -2996,6 +3045,7 @@ for (posv, potrf, potri, potrs, pstrf, elty, rtyp) in # .. Array Arguments .. # DOUBLE PRECISION A( LDA, * ), B( LDB, * ) function potrs!(uplo::AbstractChar, A::AbstractMatrix{$elty}, B::AbstractVecOrMat{$elty}) + @assert is_one_indexed(A, B) chkstride1(A, B) n = checksquare(A) chkuplo(uplo) @@ -3113,6 +3163,7 @@ for (ptsv, pttrf, elty, relty) in # .. Array Arguments .. # DOUBLE PRECISION B( LDB, * ), D( * ), E( * ) function ptsv!(D::AbstractVector{$relty}, E::AbstractVector{$elty}, B::AbstractVecOrMat{$elty}) + @assert is_one_indexed(D, E, B) chkstride1(B, D, E) n = length(D) if length(E) != n - 1 @@ -3136,6 +3187,7 @@ for (ptsv, pttrf, elty, relty) in # .. Array Arguments .. # DOUBLE PRECISION D( * ), E( * ) function pttrf!(D::AbstractVector{$relty}, E::AbstractVector{$elty}) + @assert is_one_indexed(D, E) chkstride1(D, E) n = length(D) if length(E) != n - 1 @@ -3179,6 +3231,7 @@ for (pttrs, elty, relty) in # .. Array Arguments .. # DOUBLE PRECISION B( LDB, * ), D( * ), E( * ) function pttrs!(D::AbstractVector{$relty}, E::AbstractVector{$elty}, B::AbstractVecOrMat{$elty}) + @assert is_one_indexed(D, E, B) chkstride1(B, D, E) n = length(D) if length(E) != n - 1 @@ -3211,6 +3264,7 @@ for (pttrs, elty, relty) in # DOUBLE PRECISION D( * ) # COMPLEX*16 B( LDB, * ), E( * ) function pttrs!(uplo::AbstractChar, D::AbstractVector{$relty}, E::AbstractVector{$elty}, B::AbstractVecOrMat{$elty}) + @assert is_one_indexed(D, E, B) chkstride1(B, D, E) chkuplo(uplo) n = length(D) @@ -3276,6 +3330,7 @@ for (trtri, trtrs, elty) in # DOUBLE PRECISION A( LDA, * ), B( LDB, * ) function trtrs!(uplo::AbstractChar, trans::AbstractChar, diag::AbstractChar, A::AbstractMatrix{$elty}, B::AbstractVecOrMat{$elty}) + @assert is_one_indexed(A, B) chktrans(trans) chkdiag(diag) chkstride1(A) @@ -3363,6 +3418,7 @@ for (trcon, trevc, trrfs, elty) in function trevc!(side::AbstractChar, howmny::AbstractChar, select::AbstractVector{BlasInt}, T::AbstractMatrix{$elty}, VL::AbstractMatrix{$elty} = similar(T), VR::AbstractMatrix{$elty} = similar(T)) + @assert is_one_indexed(select, T, VL, VR) # Extract if side ∉ ['L','R','B'] throw(ArgumentError("side argument must be 'L' (left eigenvectors), 'R' (right eigenvectors), or 'B' (both), got $side")) @@ -3422,6 +3478,7 @@ for (trcon, trevc, trrfs, elty) in 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))) + @assert is_one_indexed(A, B, X, Ferr, Berr) chkstride1(A, B, X, Ferr, Berr) chktrans(trans) chkuplo(uplo) @@ -3493,6 +3550,7 @@ for (trcon, trevc, trrfs, elty, relty) in function trevc!(side::AbstractChar, howmny::AbstractChar, select::AbstractVector{BlasInt}, T::AbstractMatrix{$elty}, VL::AbstractMatrix{$elty} = similar(T), VR::AbstractMatrix{$elty} = similar(T)) + @assert is_one_indexed(select, T, VL, VR) # Extract n, mm = checksquare(T), size(VL, 2) ldt, ldvl, ldvr = stride(T, 2), stride(VL, 2), stride(VR, 2) @@ -3552,6 +3610,7 @@ for (trcon, trevc, trrfs, elty, relty) in 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))) + @assert is_one_indexed(A, B, X, Ferr, Berr) chkstride1(A, B, X, Ferr, Berr) chktrans(trans) chkuplo(uplo) @@ -3626,6 +3685,7 @@ for (stev, stebz, stegr, stein, elty) in ) @eval begin function stev!(job::AbstractChar, dv::AbstractVector{$elty}, ev::AbstractVector{$elty}) + @assert is_one_indexed(dv, ev) chkstride1(dv, ev) n = length(dv) if length(ev) != n - 1 @@ -3647,6 +3707,7 @@ for (stev, stebz, stegr, stein, elty) in #* in the half-open interval (VL, VU], or the IL-th through IU-th #* eigenvalues. function stebz!(range::AbstractChar, order::AbstractChar, vl::$elty, vu::$elty, il::Integer, iu::Integer, abstol::Real, dv::AbstractVector{$elty}, ev::AbstractVector{$elty}) + @assert is_one_indexed(dv, ev) chkstride1(dv, ev) n = length(dv) if length(ev) != n - 1 @@ -3677,6 +3738,7 @@ for (stev, stebz, stegr, stein, elty) in end function stegr!(jobz::AbstractChar, range::AbstractChar, dv::AbstractVector{$elty}, ev::AbstractVector{$elty}, vl::Real, vu::Real, il::Integer, iu::Integer) + @assert is_one_indexed(dv, ev) chkstride1(dv, ev) n = length(dv) if length(ev) != n - 1 @@ -3718,6 +3780,7 @@ for (stev, stebz, stegr, stein, elty) in end function stein!(dv::AbstractVector{$elty}, ev_in::AbstractVector{$elty}, w_in::AbstractVector{$elty}, iblock_in::AbstractVector{BlasInt}, isplit_in::AbstractVector{BlasInt}) + @assert is_one_indexed(dv, ev_in, w_in, iblock_in, isplit_in) chkstride1(dv, ev_in, w_in, iblock_in, isplit_in) n = length(dv) if length(ev_in) != n - 1 @@ -3857,6 +3920,7 @@ for (syconv, sysv, sytrf, sytri, sytrs, elty) in # INTEGER IPIV( * ) # DOUBLE PRECISION A( LDA, * ), B( LDB, * ), WORK( * ) function sysv!(uplo::AbstractChar, A::AbstractMatrix{$elty}, B::AbstractVecOrMat{$elty}) + @assert is_one_indexed(A, B) chkstride1(A,B) n = checksquare(A) chkuplo(uplo) @@ -3976,6 +4040,7 @@ for (syconv, sysv, sytrf, sytri, sytrs, elty) in # DOUBLE PRECISION A( LDA, * ), B( LDB, * ) function sytrs!(uplo::AbstractChar, A::AbstractMatrix{$elty}, ipiv::AbstractVector{BlasInt}, B::AbstractVecOrMat{$elty}) + @assert is_one_indexed(A, ipiv, B) chkstride1(A,B,ipiv) n = checksquare(A) chkuplo(uplo) @@ -4007,6 +4072,7 @@ for (sysv, sytrf, sytri, sytrs, syconvf, elty) in # INTEGER IPIV( * ) # DOUBLE PRECISION A( LDA, * ), B( LDB, * ), WORK( * ) function sysv_rook!(uplo::AbstractChar, A::AbstractMatrix{$elty}, B::AbstractVecOrMat{$elty}) + @assert is_one_indexed(A, B) chkstride1(A,B) n = checksquare(A) chkuplo(uplo) @@ -4097,6 +4163,7 @@ for (sysv, sytrf, sytri, sytrs, syconvf, elty) in # DOUBLE PRECISION A( LDA, * ), B( LDB, * ) function sytrs_rook!(uplo::AbstractChar, A::AbstractMatrix{$elty}, ipiv::AbstractVector{BlasInt}, B::AbstractVecOrMat{$elty}) + @assert is_one_indexed(A, ipiv, B) chkstride1(A,B,ipiv) n = checksquare(A) chkuplo(uplo) @@ -4124,6 +4191,7 @@ for (sysv, sytrf, sytri, sytrs, syconvf, elty) in function syconvf_rook!(uplo::AbstractChar, way::AbstractChar, A::AbstractMatrix{$elty}, ipiv::AbstractVector{BlasInt}, e::AbstractVector{$elty} = Vector{$elty}(undef, length(ipiv))) + @assert is_one_indexed(A, ipiv, e) # extract n = checksquare(A) lda = max(1, stride(A, 2)) @@ -4193,6 +4261,7 @@ for (syconv, hesv, hetrf, hetri, hetrs, elty, relty) in # INTEGER IPIV( * ) # COMPLEX*16 A( LDA, * ), B( LDB, * ), WORK( * ) function hesv!(uplo::AbstractChar, A::AbstractMatrix{$elty}, B::AbstractVecOrMat{$elty}) + @assert is_one_indexed(A, B) chkstride1(A,B) n = checksquare(A) chkuplo(uplo) @@ -4310,6 +4379,7 @@ for (syconv, hesv, hetrf, hetri, hetrs, elty, relty) in # COMPLEX*16 A( LDA, * ), B( LDB, * ) function hetrs!(uplo::AbstractChar, A::AbstractMatrix{$elty}, ipiv::AbstractVector{BlasInt}, B::AbstractVecOrMat{$elty}) + @assert is_one_indexed(A, ipiv, B) chkstride1(A,B,ipiv) n = checksquare(A) if n != size(B,1) @@ -4339,6 +4409,7 @@ for (hesv, hetrf, hetri, hetrs, elty, relty) in # INTEGER IPIV( * ) # COMPLEX*16 A( LDA, * ), B( LDB, * ), WORK( * ) function hesv_rook!(uplo::AbstractChar, A::AbstractMatrix{$elty}, B::AbstractVecOrMat{$elty}) + @assert is_one_indexed(A, B) chkstride1(A,B) n = checksquare(A) chkuplo(uplo) @@ -4426,6 +4497,7 @@ for (hesv, hetrf, hetri, hetrs, elty, relty) in # COMPLEX*16 A( LDA, * ), B( LDB, * ) function hetrs_rook!(uplo::AbstractChar, A::AbstractMatrix{$elty}, ipiv::AbstractVector{BlasInt}, B::AbstractVecOrMat{$elty}) + @assert is_one_indexed(A, ipiv, B) chkstride1(A,B,ipiv) n = checksquare(A) if n != size(B,1) @@ -4456,6 +4528,7 @@ for (sysv, sytrf, sytri, sytrs, elty, relty) in # INTEGER IPIV( * ) # COMPLEX*16 A( LDA, * ), B( LDB, * ), WORK( * ) function sysv!(uplo::AbstractChar, A::AbstractMatrix{$elty}, B::AbstractVecOrMat{$elty}) + @assert is_one_indexed(A, B) chkstride1(A,B) n = checksquare(A) chkuplo(uplo) @@ -4576,6 +4649,7 @@ for (sysv, sytrf, sytri, sytrs, elty, relty) in # COMPLEX*16 A( LDA, * ), B( LDB, * ) function sytrs!(uplo::AbstractChar, A::AbstractMatrix{$elty}, ipiv::AbstractVector{BlasInt}, B::AbstractVecOrMat{$elty}) + @assert is_one_indexed(A, ipiv, B) chkstride1(A,B,ipiv) n = checksquare(A) chkuplo(uplo) @@ -4607,6 +4681,7 @@ for (sysv, sytrf, sytri, sytrs, syconvf, elty, relty) in # INTEGER IPIV( * ) # COMPLEX*16 A( LDA, * ), B( LDB, * ), WORK( * ) function sysv_rook!(uplo::AbstractChar, A::AbstractMatrix{$elty}, B::AbstractVecOrMat{$elty}) + @assert is_one_indexed(A, B) chkstride1(A,B) n = checksquare(A) chkuplo(uplo) @@ -4698,6 +4773,7 @@ for (sysv, sytrf, sytri, sytrs, syconvf, elty, relty) in # COMPLEX*16 A( LDA, * ), B( LDB, * ) function sytrs_rook!(uplo::AbstractChar, A::AbstractMatrix{$elty}, ipiv::AbstractVector{BlasInt}, B::AbstractVecOrMat{$elty}) + @assert is_one_indexed(A, ipiv, B) chkstride1(A,B,ipiv) n = checksquare(A) chkuplo(uplo) @@ -4725,6 +4801,7 @@ for (sysv, sytrf, sytri, sytrs, syconvf, elty, relty) in function syconvf_rook!(uplo::AbstractChar, way::AbstractChar, A::AbstractMatrix{$elty}, ipiv::AbstractVector{BlasInt}, e::AbstractVector{$elty} = Vector{$elty}(undef, length(ipiv))) + @assert is_one_indexed(A, ipiv, e) chkstride1(A, ipiv, e) # extract @@ -5205,6 +5282,7 @@ for (bdsqr, relty, elty) in @eval begin function bdsqr!(uplo::AbstractChar, d::AbstractVector{$relty}, e_::AbstractVector{$relty}, Vt::AbstractMatrix{$elty}, U::AbstractMatrix{$elty}, C::AbstractMatrix{$elty}) + @assert is_one_indexed(d, e_, Vt, U, C) chkstride1(d, e_, Vt, U, C) # Extract number n = length(d) @@ -5274,6 +5352,7 @@ for (bdsdc, elty) in # DOUBLE PRECISION D( * ), E( * ), Q( * ), U( LDU, * ), # $ VT( LDVT, * ), WORK( * ) function bdsdc!(uplo::AbstractChar, compq::AbstractChar, d::AbstractVector{$elty}, e_::AbstractVector{$elty}) + @assert is_one_indexed(d, e_) chkstride1(d, e_) n, ldiq, ldq, ldu, ldvt = length(d), 1, 1, 1, 1 chkuplo(uplo) @@ -5468,6 +5547,7 @@ for (orghr, elty) in # * .. Array Arguments .. # DOUBLE PRECISION A( LDA, * ), TAU( * ), WORK( * ) function orghr!(ilo::Integer, ihi::Integer, A::AbstractMatrix{$elty}, tau::AbstractVector{$elty}) + @assert is_one_indexed(A, tau) chkstride1(A, tau) n = checksquare(A) if n - length(tau) != 1 @@ -5518,6 +5598,7 @@ for (ormhr, elty) in function ormhr!(side::AbstractChar, trans::AbstractChar, ilo::Integer, ihi::Integer, A::AbstractMatrix{$elty}, tau::AbstractVector{$elty}, C::AbstractVecOrMat{$elty}) + @assert is_one_indexed(A, tau, C) chkstride1(A, tau, C) n = checksquare(A) mC, nC = size(C, 1), size(C, 2) @@ -5566,6 +5647,7 @@ for (gees, gges, elty) in # DOUBLE PRECISION A( LDA, * ), VS( LDVS, * ), WI( * ), WORK( * ), # $ WR( * ) function gees!(jobvs::AbstractChar, A::AbstractMatrix{$elty}) + @assert is_one_indexed(A) chkstride1(A) n = checksquare(A) sdim = Vector{BlasInt}(undef, 1) @@ -5659,6 +5741,7 @@ for (gees, gges, elty, relty) in # DOUBLE PRECISION RWORK( * ) # COMPLEX*16 A( LDA, * ), VS( LDVS, * ), W( * ), WORK( * ) function gees!(jobvs::AbstractChar, A::AbstractMatrix{$elty}) + @assert is_one_indexed(A) chkstride1(A) n = checksquare(A) sort = 'N' @@ -6103,6 +6186,7 @@ for (fn, elty, relty) in ((:dtrsyl_, :Float64, :Float64), @eval begin function trsyl!(transa::AbstractChar, transb::AbstractChar, A::AbstractMatrix{$elty}, B::AbstractMatrix{$elty}, C::AbstractMatrix{$elty}, isgn::Int=1) + @assert is_one_indexed(A, B, C) chkstride1(A, B, C) m, n = checksquare(A, B) lda = max(1, stride(A, 2)) diff --git a/stdlib/LinearAlgebra/src/ldlt.jl b/stdlib/LinearAlgebra/src/ldlt.jl index c97e4910014e4..ca9c1331892f2 100644 --- a/stdlib/LinearAlgebra/src/ldlt.jl +++ b/stdlib/LinearAlgebra/src/ldlt.jl @@ -92,6 +92,7 @@ end factorize(S::SymTridiagonal) = ldlt(S) function ldiv!(S::LDLt{T,M}, B::AbstractVecOrMat{T}) where {T,M<:SymTridiagonal{T}} + @assert is_one_indexed(B) n, nrhs = size(B, 1), size(B, 2) if size(S,1) != n throw(DimensionMismatch("Matrix has dimensions $(size(S)) but right hand side has first dimension $n")) diff --git a/stdlib/LinearAlgebra/src/lq.jl b/stdlib/LinearAlgebra/src/lq.jl index 57fc6b6f5d8a0..651765de608ca 100644 --- a/stdlib/LinearAlgebra/src/lq.jl +++ b/stdlib/LinearAlgebra/src/lq.jl @@ -277,6 +277,7 @@ end # With a real lhs and complex rhs with the same precision, we can reinterpret # the complex rhs as a real rhs with twice the number of columns function (\)(F::LQ{T}, B::VecOrMat{Complex{T}}) where T<:BlasReal + @assert is_one_indexed(B) c2r = reshape(copy(transpose(reinterpret(T, reshape(B, (1, length(B)))))), size(B, 1), 2*size(B, 2)) x = ldiv!(F, c2r) return reshape(copy(reinterpret(Complex{T}, copy(transpose(reshape(x, div(length(x), 2), 2))))), diff --git a/stdlib/LinearAlgebra/src/lu.jl b/stdlib/LinearAlgebra/src/lu.jl index 5426dfbd0d8ac..44c3bf46689b3 100644 --- a/stdlib/LinearAlgebra/src/lu.jl +++ b/stdlib/LinearAlgebra/src/lu.jl @@ -262,6 +262,7 @@ size(A::LU) = size(getfield(A, :factors)) size(A::LU, i) = size(getfield(A, :factors), i) function ipiv2perm(v::AbstractVector{T}, maxi::Integer) where T + @assert is_one_indexed(v) p = T[1:maxi;] @inbounds for i in 1:length(v) p[i], p[v[i]] = p[v[i]], p[i] @@ -507,6 +508,7 @@ end # See dgtts2.f function ldiv!(A::LU{T,Tridiagonal{T,V}}, B::AbstractVecOrMat) where {T,V} + @assert is_one_indexed(B) n = size(A,1) if n != size(B,1) throw(DimensionMismatch("matrix has dimensions ($n,$n) but right hand side has $(size(B,1)) rows")) @@ -538,6 +540,7 @@ function ldiv!(A::LU{T,Tridiagonal{T,V}}, B::AbstractVecOrMat) where {T,V} end function ldiv!(transA::Transpose{<:Any,<:LU{T,Tridiagonal{T,V}}}, B::AbstractVecOrMat) where {T,V} + @assert is_one_indexed(B) A = transA.parent n = size(A,1) if n != size(B,1) @@ -574,6 +577,7 @@ end # Ac_ldiv_B!(A::LU{T,Tridiagonal{T}}, B::AbstractVecOrMat) where {T<:Real} = At_ldiv_B!(A,B) function ldiv!(adjA::Adjoint{<:Any,LU{T,Tridiagonal{T,V}}}, B::AbstractVecOrMat) where {T,V} + @assert is_one_indexed(B) A = adjA.parent n = size(A,1) if n != size(B,1) diff --git a/stdlib/LinearAlgebra/src/matmul.jl b/stdlib/LinearAlgebra/src/matmul.jl index 8609ef2dc1ba3..380040d1da50a 100644 --- a/stdlib/LinearAlgebra/src/matmul.jl +++ b/stdlib/LinearAlgebra/src/matmul.jl @@ -47,7 +47,7 @@ function (*)(A::StridedMatrix{T}, x::StridedVector{S}) where {T<:BlasFloat,S} end function (*)(A::AbstractMatrix{T}, x::AbstractVector{S}) where {T,S} TS = promote_op(matprod, T, S) - mul!(similar(x,TS,size(A,1)),A,x) + mul!(similar(x,TS,axes(A,1)),A,x) end # these will throw a DimensionMismatch unless B has 1 row (or 1 col for transposed case): @@ -492,6 +492,7 @@ end # strides != 1 cases function generic_matvecmul!(C::AbstractVector{R}, tA, A::AbstractVecOrMat, B::AbstractVector) where R + @assert is_one_indexed(C, A, B) mB = length(B) mA, nA = lapack_size(tA, A) if mB != nA @@ -578,6 +579,7 @@ end generic_matmatmul!(C::AbstractVecOrMat, tA, tB, A::AbstractVecOrMat, B::AbstractVecOrMat) = _generic_matmatmul!(C, tA, tB, A, B) function _generic_matmatmul!(C::AbstractVecOrMat{R}, tA, tB, A::AbstractVecOrMat{T}, B::AbstractVecOrMat{S}) where {T,S,R} + @assert is_one_indexed(C, A, B) mA, nA = lapack_size(tA, A) mB, nB = lapack_size(tB, B) if mB != nA @@ -751,6 +753,7 @@ function matmul2x2(tA, tB, A::AbstractMatrix{T}, B::AbstractMatrix{S}) where {T, end function matmul2x2!(C::AbstractMatrix, tA, tB, A::AbstractMatrix, B::AbstractMatrix) + @assert is_one_indexed(C, A, B) if !(size(A) == size(B) == size(C) == (2,2)) throw(DimensionMismatch("A has size $(size(A)), B has size $(size(B)), C has size $(size(C))")) end @@ -792,6 +795,7 @@ function matmul3x3(tA, tB, A::AbstractMatrix{T}, B::AbstractMatrix{S}) where {T, end function matmul3x3!(C::AbstractMatrix, tA, tB, A::AbstractMatrix, B::AbstractMatrix) + @assert is_one_indexed(C, A, B) if !(size(A) == size(B) == size(C) == (3,3)) throw(DimensionMismatch("A has size $(size(A)), B has size $(size(B)), C has size $(size(C))")) end diff --git a/stdlib/LinearAlgebra/src/qr.jl b/stdlib/LinearAlgebra/src/qr.jl index 602c143230cfb..793a66033293f 100644 --- a/stdlib/LinearAlgebra/src/qr.jl +++ b/stdlib/LinearAlgebra/src/qr.jl @@ -156,6 +156,7 @@ Base.iterate(S::QRPivoted, ::Val{:p}) = (S.p, Val(:done)) Base.iterate(S::QRPivoted, ::Val{:done}) = nothing function qrfactUnblocked!(A::AbstractMatrix{T}) where {T} + @assert is_one_indexed(A) m, n = size(A) τ = zeros(T, min(m,n)) for k = 1:min(m - 1 + !(T<:Real), n) @@ -332,17 +333,22 @@ true compactly rather as two separate dense matrices. """ function qr(A::AbstractMatrix{T}, arg) where T + @assert is_one_indexed(A) AA = similar(A, _qreltype(T), size(A)) copyto!(AA, A) return qr!(AA, arg) end function qr(A::AbstractMatrix{T}) where T + @assert is_one_indexed(A) AA = similar(A, _qreltype(T), size(A)) copyto!(AA, A) return qr!(AA) end qr(x::Number) = qr(fill(x,1,1)) -qr(v::AbstractVector) = qr(reshape(v, (length(v), 1))) +function qr(v::AbstractVector) + @assert is_one_indexed(v) + qr(reshape(v, (length(v), 1))) +end # Conversions QR{T}(A::QR) where {T} = QR(convert(AbstractMatrix{T}, A.factors), convert(Vector{T}, A.τ)) @@ -479,6 +485,7 @@ lmul!(A::QRCompactWYQ{T,S}, B::StridedVecOrMat{T}) where {T<:BlasFloat, S<:Strid lmul!(A::QRPackedQ{T,S}, B::StridedVecOrMat{T}) where {T<:BlasFloat, S<:StridedMatrix} = LAPACK.ormqr!('L','N',A.factors,A.τ,B) function lmul!(A::QRPackedQ, B::AbstractVecOrMat) + @assert is_one_indexed(B) mA, nA = size(A.factors) mB, nB = size(B,1), size(B,2) if mA != mB @@ -538,6 +545,7 @@ lmul!(adjA::Adjoint{<:Any,<:QRPackedQ{T,S}}, B::StridedVecOrMat{T}) where {T<:Bl lmul!(adjA::Adjoint{<:Any,<:QRPackedQ{T,S}}, B::StridedVecOrMat{T}) where {T<:BlasComplex,S<:StridedMatrix} = (A = adjA.parent; LAPACK.ormqr!('L','C',A.factors,A.τ,B)) function lmul!(adjA::Adjoint{<:Any,<:QRPackedQ}, B::AbstractVecOrMat) + @assert is_one_indexed(B) A = adjA.parent mA, nA = size(A.factors) mB, nB = size(B,1), size(B,2) @@ -801,6 +809,7 @@ _zeros(::Type{T}, b::AbstractVector, n::Integer) where {T} = zeros(T, max(length _zeros(::Type{T}, B::AbstractMatrix, n::Integer) where {T} = zeros(T, max(size(B, 1), n), size(B, 2)) function (\)(A::Union{QR{TA},QRCompactWY{TA},QRPivoted{TA}}, B::AbstractVecOrMat{TB}) where {TA,TB} + @assert is_one_indexed(B) S = promote_type(TA,TB) m, n = size(A) m == size(B,1) || throw(DimensionMismatch("left hand side has $m rows, but right hand side has $(size(B,1)) rows")) @@ -821,6 +830,7 @@ _ret_size(A::Factorization, b::AbstractVector) = (max(size(A, 2), length(b)),) _ret_size(A::Factorization, B::AbstractMatrix) = (max(size(A, 2), size(B, 1)), size(B, 2)) function (\)(A::Union{QR{T},QRCompactWY{T},QRPivoted{T}}, BIn::VecOrMat{Complex{T}}) where T<:BlasReal + @assert is_one_indexed(BIn) m, n = size(A) m == size(BIn, 1) || throw(DimensionMismatch("left hand side has $m rows, but right hand side has $(size(BIn,1)) rows")) diff --git a/stdlib/LinearAlgebra/src/triangular.jl b/stdlib/LinearAlgebra/src/triangular.jl index 029501cc49bd2..5875c3867a2b9 100644 --- a/stdlib/LinearAlgebra/src/triangular.jl +++ b/stdlib/LinearAlgebra/src/triangular.jl @@ -1138,6 +1138,7 @@ end # replacing repeated references to A.data[j,j] with [Ajj = A.data[j,j] and references to Ajj] # does not significantly impact performance as of Dec 2015 function naivesub!(A::UpperTriangular, b::AbstractVector, x::AbstractVector = b) + @assert is_one_indexed(A, b, x) n = size(A, 2) if !(n == length(b) == length(x)) throw(DimensionMismatch("second dimension of left hand side A, $n, length of output x, $(length(x)), and length of right hand side b, $(length(b)), must be equal")) @@ -1152,6 +1153,7 @@ function naivesub!(A::UpperTriangular, b::AbstractVector, x::AbstractVector = b) x end function naivesub!(A::UnitUpperTriangular, b::AbstractVector, x::AbstractVector = b) + @assert is_one_indexed(A, b, x) n = size(A, 2) if !(n == length(b) == length(x)) throw(DimensionMismatch("second dimension of left hand side A, $n, length of output x, $(length(x)), and length of right hand side b, $(length(b)), must be equal")) @@ -1165,6 +1167,7 @@ function naivesub!(A::UnitUpperTriangular, b::AbstractVector, x::AbstractVector x end function naivesub!(A::LowerTriangular, b::AbstractVector, x::AbstractVector = b) + @assert is_one_indexed(A, b, x) n = size(A, 2) if !(n == length(b) == length(x)) throw(DimensionMismatch("second dimension of left hand side A, $n, length of output x, $(length(x)), and length of right hand side b, $(length(b)), must be equal")) @@ -1179,6 +1182,7 @@ function naivesub!(A::LowerTriangular, b::AbstractVector, x::AbstractVector = b) x end function naivesub!(A::UnitLowerTriangular, b::AbstractVector, x::AbstractVector = b) + @assert is_one_indexed(A, b, x) n = size(A, 2) if !(n == length(b) == length(x)) throw(DimensionMismatch("second dimension of left hand side A, $n, length of output x, $(length(x)), and length of right hand side b, $(length(b)), must be equal")) @@ -1194,6 +1198,7 @@ end # in the following transpose and conjugate transpose naive substitution variants, # accumulating in z rather than b[j] significantly improves performance as of Dec 2015 function ldiv!(transA::Transpose{<:Any,<:LowerTriangular}, b::AbstractVector, x::AbstractVector) + @assert is_one_indexed(transA, b, x) A = transA.parent n = size(A, 1) if !(n == length(b) == length(x)) @@ -1212,6 +1217,7 @@ end ldiv!(transA::Transpose{<:Any,<:LowerTriangular}, b::AbstractVector) = ldiv!(transA, b, b) function ldiv!(transA::Transpose{<:Any,<:UnitLowerTriangular}, b::AbstractVector, x::AbstractVector) + @assert is_one_indexed(transA, b, x) A = transA.parent n = size(A, 1) if !(n == length(b) == length(x)) @@ -1229,6 +1235,7 @@ end ldiv!(transA::Transpose{<:Any,<:UnitLowerTriangular}, b::AbstractVector) = ldiv!(transA, b, b) function ldiv!(transA::Transpose{<:Any,<:UpperTriangular}, b::AbstractVector, x::AbstractVector) + @assert is_one_indexed(transA, b, x) A = transA.parent n = size(A, 1) if !(n == length(b) == length(x)) @@ -1247,6 +1254,7 @@ end ldiv!(transA::Transpose{<:Any,<:UpperTriangular}, b::AbstractVector) = ldiv!(transA, b, b) function ldiv!(transA::Transpose{<:Any,<:UnitUpperTriangular}, b::AbstractVector, x::AbstractVector) + @assert is_one_indexed(transA, b, x) A = transA.parent n = size(A, 1) if !(n == length(b) == length(x)) @@ -1264,6 +1272,7 @@ end ldiv!(transA::Transpose{<:Any,<:UnitUpperTriangular}, b::AbstractVector) = ldiv!(transA, b, b) function ldiv!(adjA::Adjoint{<:Any,<:LowerTriangular}, b::AbstractVector, x::AbstractVector) + @assert is_one_indexed(adjA, b, x) A = adjA.parent n = size(A, 1) if !(n == length(b) == length(x)) @@ -1282,6 +1291,7 @@ end ldiv!(adjA::Adjoint{<:Any,<:LowerTriangular}, b::AbstractVector) = ldiv!(adjA, b, b) function ldiv!(adjA::Adjoint{<:Any,<:UnitLowerTriangular}, b::AbstractVector, x::AbstractVector) + @assert is_one_indexed(adjA, b, x) A = adjA.parent n = size(A, 1) if !(n == length(b) == length(x)) @@ -1299,6 +1309,7 @@ end ldiv!(adjA::Adjoint{<:Any,<:UnitLowerTriangular}, b::AbstractVector) = ldiv!(adjA, b, b) function ldiv!(adjA::Adjoint{<:Any,<:UpperTriangular}, b::AbstractVector, x::AbstractVector) + @assert is_one_indexed(adjA, b, x) A = adjA.parent n = size(A, 1) if !(n == length(b) == length(x)) @@ -1317,6 +1328,7 @@ end ldiv!(adjA::Adjoint{<:Any,<:UpperTriangular}, b::AbstractVector) = ldiv!(adjA, b, b) function ldiv!(adjA::Adjoint{<:Any,<:UnitUpperTriangular}, b::AbstractVector, x::AbstractVector) + @assert is_one_indexed(adjA, b, x) A = adjA.parent n = size(A, 1) if !(n == length(b) == length(x)) @@ -1783,12 +1795,14 @@ for mat in (:AbstractVector, :AbstractMatrix) ### Multiplication with triangle to the left and hence rhs cannot be transposed. @eval begin function *(A::AbstractTriangular, B::$mat) + @assert is_one_indexed(B) TAB = typeof(zero(eltype(A))*zero(eltype(B)) + zero(eltype(A))*zero(eltype(B))) BB = similar(B, TAB, size(B)) copyto!(BB, B) lmul!(convert(AbstractArray{TAB}, A), BB) end function *(adjA::Adjoint{<:Any,<:AbstractTriangular}, B::$mat) + @assert is_one_indexed(B) A = adjA.parent TAB = typeof(zero(eltype(A))*zero(eltype(B)) + zero(eltype(A))*zero(eltype(B))) BB = similar(B, TAB, size(B)) @@ -1796,6 +1810,7 @@ for mat in (:AbstractVector, :AbstractMatrix) lmul!(adjoint(convert(AbstractArray{TAB}, A)), BB) end function *(transA::Transpose{<:Any,<:AbstractTriangular}, B::$mat) + @assert is_one_indexed(B) A = transA.parent TAB = typeof(zero(eltype(A))*zero(eltype(B)) + zero(eltype(A))*zero(eltype(B))) BB = similar(B, TAB, size(B)) @@ -1806,12 +1821,14 @@ for mat in (:AbstractVector, :AbstractMatrix) ### Left division with triangle to the left hence rhs cannot be transposed. No quotients. @eval begin function \(A::Union{UnitUpperTriangular,UnitLowerTriangular}, B::$mat) + @assert is_one_indexed(B) TAB = typeof(zero(eltype(A))*zero(eltype(B)) + zero(eltype(A))*zero(eltype(B))) BB = similar(B, TAB, size(B)) copyto!(BB, B) ldiv!(convert(AbstractArray{TAB}, A), BB) end function \(adjA::Adjoint{<:Any,<:Union{UnitUpperTriangular,UnitLowerTriangular}}, B::$mat) + @assert is_one_indexed(B) A = adjA.parent TAB = typeof(zero(eltype(A))*zero(eltype(B)) + zero(eltype(A))*zero(eltype(B))) BB = similar(B, TAB, size(B)) @@ -1819,6 +1836,7 @@ for mat in (:AbstractVector, :AbstractMatrix) ldiv!(adjoint(convert(AbstractArray{TAB}, A)), BB) end function \(transA::Transpose{<:Any,<:Union{UnitUpperTriangular,UnitLowerTriangular}}, B::$mat) + @assert is_one_indexed(B) A = transA.parent TAB = typeof(zero(eltype(A))*zero(eltype(B)) + zero(eltype(A))*zero(eltype(B))) BB = similar(B, TAB, size(B)) @@ -1829,12 +1847,14 @@ for mat in (:AbstractVector, :AbstractMatrix) ### Left division with triangle to the left hence rhs cannot be transposed. Quotients. @eval begin function \(A::Union{UpperTriangular,LowerTriangular}, B::$mat) + @assert is_one_indexed(B) TAB = typeof((zero(eltype(A))*zero(eltype(B)) + zero(eltype(A))*zero(eltype(B)))/one(eltype(A))) BB = similar(B, TAB, size(B)) copyto!(BB, B) ldiv!(convert(AbstractArray{TAB}, A), BB) end function \(adjA::Adjoint{<:Any,<:Union{UpperTriangular,LowerTriangular}}, B::$mat) + @assert is_one_indexed(B) A = adjA.parent TAB = typeof((zero(eltype(A))*zero(eltype(B)) + zero(eltype(A))*zero(eltype(B)))/one(eltype(A))) BB = similar(B, TAB, size(B)) @@ -1842,6 +1862,7 @@ for mat in (:AbstractVector, :AbstractMatrix) ldiv!(adjoint(convert(AbstractArray{TAB}, A)), BB) end function \(transA::Transpose{<:Any,<:Union{UpperTriangular,LowerTriangular}}, B::$mat) + @assert is_one_indexed(B) A = transA.parent TAB = typeof((zero(eltype(A))*zero(eltype(B)) + zero(eltype(A))*zero(eltype(B)))/one(eltype(A))) BB = similar(B, TAB, size(B)) @@ -1851,6 +1872,7 @@ for mat in (:AbstractVector, :AbstractMatrix) end ### Right division with triangle to the right hence lhs cannot be transposed. No quotients. @eval begin + @assert is_one_indexed(A) function /(A::$mat, B::Union{UnitUpperTriangular, UnitLowerTriangular}) TAB = typeof(zero(eltype(A))*zero(eltype(B)) + zero(eltype(A))*zero(eltype(B))) AA = similar(A, TAB, size(A)) @@ -1858,6 +1880,7 @@ for mat in (:AbstractVector, :AbstractMatrix) rdiv!(AA, convert(AbstractArray{TAB}, B)) end function /(A::$mat, adjB::Adjoint{<:Any,<:Union{UnitUpperTriangular, UnitLowerTriangular}}) + @assert is_one_indexed(A) B = adjB.parent TAB = typeof(zero(eltype(A))*zero(eltype(B)) + zero(eltype(A))*zero(eltype(B))) AA = similar(A, TAB, size(A)) @@ -1865,6 +1888,7 @@ for mat in (:AbstractVector, :AbstractMatrix) rdiv!(AA, adjoint(convert(AbstractArray{TAB}, B))) end function /(A::$mat, transB::Transpose{<:Any,<:Union{UnitUpperTriangular, UnitLowerTriangular}}) + @assert is_one_indexed(A) B = transB.parent TAB = typeof(zero(eltype(A))*zero(eltype(B)) + zero(eltype(A))*zero(eltype(B))) AA = similar(A, TAB, size(A)) @@ -1875,12 +1899,14 @@ for mat in (:AbstractVector, :AbstractMatrix) ### Right division with triangle to the right hence lhs cannot be transposed. Quotients. @eval begin function /(A::$mat, B::Union{UpperTriangular,LowerTriangular}) + @assert is_one_indexed(A) TAB = typeof((zero(eltype(A))*zero(eltype(B)) + zero(eltype(A))*zero(eltype(B)))/one(eltype(A))) AA = similar(A, TAB, size(A)) copyto!(AA, A) rdiv!(AA, convert(AbstractArray{TAB}, B)) end function /(A::$mat, adjB::Adjoint{<:Any,<:Union{UpperTriangular,LowerTriangular}}) + @assert is_one_indexed(A) B = adjB.parent TAB = typeof((zero(eltype(A))*zero(eltype(B)) + zero(eltype(A))*zero(eltype(B)))/one(eltype(A))) AA = similar(A, TAB, size(A)) @@ -1888,6 +1914,7 @@ for mat in (:AbstractVector, :AbstractMatrix) rdiv!(AA, adjoint(convert(AbstractArray{TAB}, B))) end function /(A::$mat, transB::Transpose{<:Any,<:Union{UpperTriangular,LowerTriangular}}) + @assert is_one_indexed(A) B = transB.parent TAB = typeof((zero(eltype(A))*zero(eltype(B)) + zero(eltype(A))*zero(eltype(B)))/one(eltype(A))) AA = similar(A, TAB, size(A)) @@ -1899,12 +1926,14 @@ end ### Multiplication with triangle to the right and hence lhs cannot be transposed. # Only for AbstractMatrix, hence outside the above loop. function *(A::AbstractMatrix, B::AbstractTriangular) + @assert is_one_indexed(A) TAB = typeof(zero(eltype(A))*zero(eltype(B)) + zero(eltype(A))*zero(eltype(B))) AA = similar(A, TAB, size(A)) copyto!(AA, A) rmul!(AA, convert(AbstractArray{TAB}, B)) end function *(A::AbstractMatrix, adjB::Adjoint{<:Any,<:AbstractTriangular}) + @assert is_one_indexed(A) B = adjB.parent TAB = typeof(zero(eltype(A))*zero(eltype(B)) + zero(eltype(A))*zero(eltype(B))) AA = similar(A, TAB, size(A)) @@ -1912,6 +1941,7 @@ function *(A::AbstractMatrix, adjB::Adjoint{<:Any,<:AbstractTriangular}) rmul!(AA, adjoint(convert(AbstractArray{TAB}, B))) end function *(A::AbstractMatrix, transB::Transpose{<:Any,<:AbstractTriangular}) + @assert is_one_indexed(A) B = transB.parent TAB = typeof(zero(eltype(A))*zero(eltype(B)) + zero(eltype(A))*zero(eltype(B))) AA = similar(A, TAB, size(A)) @@ -2227,6 +2257,7 @@ end # Repeatedly compute the square roots of A so that in the end its # eigenvalues are close enough to the positive real line function invsquaring(A0::UpperTriangular, theta) + @assert is_one_indexed(theta) # assumes theta is in ascending order maxsqrt = 100 tmax = size(theta, 1) diff --git a/stdlib/LinearAlgebra/src/tridiag.jl b/stdlib/LinearAlgebra/src/tridiag.jl index a5e7a053184af..d5e853b752eaa 100644 --- a/stdlib/LinearAlgebra/src/tridiag.jl +++ b/stdlib/LinearAlgebra/src/tridiag.jl @@ -314,6 +314,7 @@ end # Linear Algebra and its Applications 212-213 (1994), pp.413-414 # doi:10.1016/0024-3795(94)90414-6 function inv_usmani(a::V, b::V, c::V) where {T,V<:AbstractVector{T}} + @assert is_one_indexed(a, b, c) n = length(b) θ = zeros(T, n+1) #principal minors of A θ[1] = 1 @@ -344,6 +345,7 @@ end #Implements the determinant using principal minors #Inputs and reference are as above for inv_usmani() function det_usmani(a::V, b::V, c::V) where {T,V<:AbstractVector{T}} + @assert is_one_indexed(a, b, c) n = length(b) θa = one(T) if n == 0 @@ -391,6 +393,7 @@ struct Tridiagonal{T,V<:AbstractVector{T}} <: AbstractMatrix{T} du::V # sup-diagonal du2::V # supsup-diagonal for pivoting in LU function Tridiagonal{T}(dl::V, d::V, du::V) where {T,V<:AbstractVector{T}} + @assert is_one_indexed(dl, d, du) n = length(d) if (length(dl) != n-1 || length(du) != n-1) throw(ArgumentError(string("cannot construct Tridiagonal from incompatible ", @@ -401,6 +404,8 @@ struct Tridiagonal{T,V<:AbstractVector{T}} <: AbstractMatrix{T} end # constructor used in lu! function Tridiagonal{T,V}(dl::V, d::V, du::V, du2::V) where {T,V<:AbstractVector{T}} + @assert is_one_indexed(dl, d, du, du2) + # length checks? new{T,V}(dl, d, du, du2) end end diff --git a/stdlib/LinearAlgebra/src/uniformscaling.jl b/stdlib/LinearAlgebra/src/uniformscaling.jl index 67d2456a5c316..e5383bfc88dd3 100644 --- a/stdlib/LinearAlgebra/src/uniformscaling.jl +++ b/stdlib/LinearAlgebra/src/uniformscaling.jl @@ -30,6 +30,8 @@ struct UniformScaling{T<:Number} λ::T end +Base.is_one_indexed(::UniformScaling) = true + """ I @@ -172,6 +174,7 @@ Broadcast.broadcasted(::typeof(/), J::UniformScaling,x::Number) = UniformScaling ## equality comparison with UniformScaling ==(J::UniformScaling, A::AbstractMatrix) = A == J function ==(A::AbstractMatrix, J::UniformScaling) + @assert is_one_indexed(A) size(A, 1) == size(A, 2) || return false iszero(J.λ) && return iszero(A) isone(J.λ) && return isone(A) @@ -204,6 +207,7 @@ end isapprox(A::AbstractMatrix, J::UniformScaling; kwargs...) = isapprox(J, A; kwargs...) function copyto!(A::AbstractMatrix, J::UniformScaling) + @assert is_one_indexed(A) size(A,1)==size(A,2) || throw(DimensionMismatch("a UniformScaling can only be copied to a square matrix")) fill!(A, 0) λ = J.λ @@ -240,6 +244,7 @@ for (f,dim,name) in ((:hcat,1,"rows"), (:vcat,2,"cols")) n = 0 for a in A if !isa(a, UniformScaling) + @assert is_one_indexed(a) na = size(a,$dim) n > 0 && n != na && throw(DimensionMismatch(string("number of ", $name, @@ -255,6 +260,7 @@ end function hvcat(rows::Tuple{Vararg{Int}}, A::Union{AbstractVecOrMat,UniformScaling}...) + @assert is_one_indexed(A...) nr = length(rows) sum(rows) == length(A) || throw(ArgumentError("mismatch between row sizes and number of arguments")) n = zeros(Int, length(A)) diff --git a/stdlib/Random/src/misc.jl b/stdlib/Random/src/misc.jl index 74d161d32b498..ec24839093312 100644 --- a/stdlib/Random/src/misc.jl +++ b/stdlib/Random/src/misc.jl @@ -85,6 +85,7 @@ end # (Note that this is different from the problem of finding a random # size-m subset of A where m is fixed!) function randsubseq!(r::AbstractRNG, S::AbstractArray, A::AbstractArray, p::Real) + @assert is_one_indexed(S, A) 0 <= p <= 1 || throw(ArgumentError("probability $p not in [0,1]")) n = length(A) p == 1 && return copyto!(resize!(S, n), A) @@ -179,6 +180,7 @@ julia> shuffle!(rng, Vector(1:16)) ``` """ function shuffle!(r::AbstractRNG, a::AbstractArray) + @assert is_one_indexed(a) n = length(a) @assert n <= Int64(2)^52 mask = nextpow2(n) - 1 diff --git a/stdlib/Serialization/src/Serialization.jl b/stdlib/Serialization/src/Serialization.jl index e60916611fbf9..fe3b57ab36926 100644 --- a/stdlib/Serialization/src/Serialization.jl +++ b/stdlib/Serialization/src/Serialization.jl @@ -230,6 +230,7 @@ function serialize(s::AbstractSerializer, x::Symbol) end function serialize_array_data(s::IO, a) + @assert is_one_indexed(a) isempty(a) && return 0 if eltype(a) === Bool last = a[1] diff --git a/stdlib/SparseArrays/src/linalg.jl b/stdlib/SparseArrays/src/linalg.jl index cc765ea42a863..55a90ea46a4e7 100644 --- a/stdlib/SparseArrays/src/linalg.jl +++ b/stdlib/SparseArrays/src/linalg.jl @@ -237,6 +237,7 @@ end ## solvers function fwdTriSolve!(A::SparseMatrixCSCUnion, B::AbstractVecOrMat) # forward substitution for CSC matrices + @assert is_one_indexed(A, B) nrowB, ncolB = size(B, 1), size(B, 2) ncol = LinearAlgebra.checksquare(A) if nrowB != ncol @@ -282,6 +283,7 @@ end function bwdTriSolve!(A::SparseMatrixCSCUnion, B::AbstractVecOrMat) # backward substitution for CSC matrices + @assert is_one_indexed(A, B) nrowB, ncolB = size(B, 1), size(B, 2) ncol = LinearAlgebra.checksquare(A) if nrowB != ncol @@ -920,6 +922,7 @@ function lmul!(b::Number, A::SparseMatrixCSC) end function \(A::SparseMatrixCSC, B::AbstractVecOrMat) + @assert is_one_indexed(A, B) m, n = size(A) if m == n if istril(A) @@ -943,6 +946,7 @@ for (xformtype, xformop) in ((:Adjoint, :adjoint), (:Transpose, :transpose)) @eval begin function \(xformA::($xformtype){<:Any,<:SparseMatrixCSC}, B::AbstractVecOrMat) A = xformA.parent + @assert is_one_indexed(A, B) m, n = size(A) if m == n if istril(A) diff --git a/stdlib/SparseArrays/src/sparsematrix.jl b/stdlib/SparseArrays/src/sparsematrix.jl index a02359d981a3b..c3fbd0a659806 100644 --- a/stdlib/SparseArrays/src/sparsematrix.jl +++ b/stdlib/SparseArrays/src/sparsematrix.jl @@ -371,6 +371,7 @@ SparseMatrixCSC(M::Matrix) = sparse(M) SparseMatrixCSC(M::AbstractMatrix{Tv}) where {Tv} = SparseMatrixCSC{Tv,Int}(M) SparseMatrixCSC{Tv}(M::AbstractMatrix{Tv}) where {Tv} = SparseMatrixCSC{Tv,Int}(M) function SparseMatrixCSC{Tv,Ti}(M::AbstractMatrix) where {Tv,Ti} + @assert is_one_indexed(M) I = findall(x -> x != 0, M) eltypeTiI = Ti[i[1] for i in I] eltypeTiJ = Ti[i[2] for i in I] @@ -449,6 +450,7 @@ sparse_IJ_sorted!(I,J,V::AbstractVector{Bool},m,n) = sparse_IJ_sorted!(I,J,V,m,n function sparse_IJ_sorted!(I::AbstractVector{Ti}, J::AbstractVector{Ti}, V::AbstractVector, m::Integer, n::Integer, combine::Function) where Ti<:Integer + @assert is_one_indexed(I, J, V) m = m < 0 ? 0 : m n = n < 0 ? 0 : n if isempty(V); return spzeros(eltype(V),Ti,m,n); end @@ -520,6 +522,7 @@ julia> sparse(Is, Js, Vs) ``` """ function sparse(I::AbstractVector{Ti}, J::AbstractVector{Ti}, V::AbstractVector{Tv}, m::Integer, n::Integer, combine) where {Tv,Ti<:Integer} + @assert is_one_indexed(I, J, V) coolen = length(I) if length(J) != coolen || length(V) != coolen throw(ArgumentError(string("the first three arguments' lengths must match, ", @@ -613,6 +616,7 @@ function sparse!(I::AbstractVector{Ti}, J::AbstractVector{Ti}, csrrowptr::Vector{Ti}, csrcolval::Vector{Ti}, csrnzval::Vector{Tv}, csccolptr::Vector{Ti}, cscrowval::Vector{Ti}, cscnzval::Vector{Tv}) where {Tv,Ti<:Integer} + @assert is_one_indexed(I, J, V) # Compute the CSR form's row counts and store them shifted forward by one in csrrowptr fill!(csrrowptr, Ti(0)) coolen = length(I) @@ -943,7 +947,8 @@ column-permutation argument `q`. """ function _checkargs_sourcecompatperms_permute!(A::SparseMatrixCSC, p::AbstractVector{<:Integer}, q::AbstractVector{<:Integer}) - if length(q) != A.n + @assert is_one_indexed(p, q) + if length(q) != A.n throw(DimensionMismatch(string("the length of column-permutation argument `q`, ", "`length(q) (= $(length(q)))`, must match source argument `A`'s column ", "count, `A.n (= $(A.n))`"))) @@ -968,6 +973,7 @@ function _checkargs_permutationsvalid_permute!( end function _ispermutationvalid_permute!(perm::AbstractVector{<:Integer}, checkspace::Vector{<:Integer}) + @assert is_one_indexed(perm) n = length(perm) checkspace[1:n] .= 0 for k in perm @@ -1665,6 +1671,7 @@ end # General mapreducedim function _mapreducerows!(f, op, R::AbstractArray, A::SparseMatrixCSC{T}) where T + @assert is_one_indexed(A, R) colptr = A.colptr rowval = A.rowval nzval = A.nzval @@ -1680,6 +1687,7 @@ function _mapreducerows!(f, op, R::AbstractArray, A::SparseMatrixCSC{T}) where T end function _mapreducecols!(f, op, R::AbstractArray, A::SparseMatrixCSC{Tv,Ti}) where {Tv,Ti} + @assert is_one_indexed(A, R) colptr = A.colptr rowval = A.rowval nzval = A.nzval @@ -1699,6 +1707,7 @@ function _mapreducecols!(f, op, R::AbstractArray, A::SparseMatrixCSC{Tv,Ti}) whe end function Base._mapreducedim!(f, op, R::AbstractArray, A::SparseMatrixCSC{T}) where T + @assert is_one_indexed(A, R) lsiz = Base.check_reducedims(R,A) isempty(A) && return R @@ -1749,6 +1758,7 @@ end # Specialized mapreducedim for + cols to avoid allocating a # temporary array when f(0) == 0 function _mapreducecols!(f, op::typeof(+), R::AbstractArray, A::SparseMatrixCSC{Tv,Ti}) where {Tv,Ti} + @assert is_one_indexed(A, R) nzval = A.nzval m, n = size(A) if length(nzval) == m*n @@ -1814,6 +1824,7 @@ function _findz(A::SparseMatrixCSC{Tv,Ti}, rows=1:A.m, cols=1:A.n) where {Tv,Ti} end function _findr(op, A, region, Tv) + @assert is_one_indexed(A) Ti = eltype(keys(A)) i1 = first(keys(A)) N = nnz(A) @@ -1920,6 +1931,7 @@ getindex(A::SparseMatrixCSC, i, ::Colon) = getindex(A, i, 1:size(A, 2)) getindex(A::SparseMatrixCSC, ::Colon, i) = getindex(A, 1:size(A, 1), i) function getindex_cols(A::SparseMatrixCSC{Tv,Ti}, J::AbstractVector) where {Tv,Ti} + @assert is_one_indexed(A, J) # for indexing whole columns (m, n) = size(A) nJ = length(J) @@ -1956,6 +1968,7 @@ getindex_traverse_col(::AbstractUnitRange, lo::Int, hi::Int) = lo:hi getindex_traverse_col(I::StepRange, lo::Int, hi::Int) = step(I) > 0 ? (lo:1:hi) : (hi:-1:lo) function getindex(A::SparseMatrixCSC{Tv,Ti}, I::AbstractRange, J::AbstractVector) where {Tv,Ti<:Integer} + @assert is_one_indexed(A, I, J) # Ranges for indexing rows (m, n) = size(A) # whole columns: @@ -2003,6 +2016,7 @@ function getindex(A::SparseMatrixCSC{Tv,Ti}, I::AbstractRange, J::AbstractVector end function getindex_I_sorted(A::SparseMatrixCSC{Tv,Ti}, I::AbstractVector, J::AbstractVector) where {Tv,Ti} + @assert is_one_indexed(A, I, J) # Sorted vectors for indexing rows. # Similar to getindex_general but without the transpose trick. (m, n) = size(A) @@ -2023,6 +2037,7 @@ function getindex_I_sorted(A::SparseMatrixCSC{Tv,Ti}, I::AbstractVector, J::Abst end function getindex_I_sorted_bsearch_A(A::SparseMatrixCSC{Tv,Ti}, I::AbstractVector, J::AbstractVector) where {Tv,Ti} + @assert is_one_indexed(A, I, J) nI = length(I) nJ = length(J) @@ -2082,6 +2097,7 @@ function getindex_I_sorted_bsearch_A(A::SparseMatrixCSC{Tv,Ti}, I::AbstractVecto end function getindex_I_sorted_linear(A::SparseMatrixCSC{Tv,Ti}, I::AbstractVector, J::AbstractVector) where {Tv,Ti} + @assert is_one_indexed(A, I, J) nI = length(I) nJ = length(J) @@ -2141,6 +2157,7 @@ function getindex_I_sorted_linear(A::SparseMatrixCSC{Tv,Ti}, I::AbstractVector, end function getindex_I_sorted_bsearch_I(A::SparseMatrixCSC{Tv,Ti}, I::AbstractVector, J::AbstractVector) where {Tv,Ti} + @assert is_one_indexed(A, I, J) nI = length(I) nJ = length(J) @@ -2242,6 +2259,7 @@ function permute_rows!(S::SparseMatrixCSC{Tv,Ti}, pI::Vector{Int}) where {Tv,Ti} end function getindex_general(A::SparseMatrixCSC, I::AbstractVector, J::AbstractVector) + @assert is_one_indexed(A, I, J) pI = sortperm(I) @inbounds Is = I[pI] permute_rows!(getindex_I_sorted(A, Is, J), pI) @@ -2249,6 +2267,7 @@ end # the general case: function getindex(A::SparseMatrixCSC{Tv,Ti}, I::AbstractVector, J::AbstractVector) where {Tv,Ti} + @assert is_one_indexed(A, I, J) (m, n) = size(A) if !isempty(J) @@ -2273,6 +2292,7 @@ function getindex(A::SparseMatrixCSC{Tv,Ti}, I::AbstractVector, J::AbstractVecto end function getindex(A::SparseMatrixCSC{Tv,Ti}, I::AbstractArray) where {Tv,Ti} + @assert is_one_indexed(A, I) szA = size(A) nA = szA[1]*szA[2] colptrA = A.colptr @@ -2375,6 +2395,7 @@ j in J, assigns zero to A[i,j] if A[i,j] is a presently-stored entry, and otherw """ function _spsetz_setindex!(A::SparseMatrixCSC, I::Union{Integer, AbstractVector{<:Integer}}, J::Union{Integer, AbstractVector{<:Integer}}) + @assert is_one_indexed(A, I, J) lengthI = length(I) for j in J coljAfirstk = A.colptr[j] @@ -2411,6 +2432,7 @@ assigns x to A[i,j] if A[i,j] is not presently stored. """ function _spsetnz_setindex!(A::SparseMatrixCSC{Tv}, x::Tv, I::Union{Integer, AbstractVector{<:Integer}}, J::Union{Integer, AbstractVector{<:Integer}}) where Tv + @assert is_one_indexed(A, I, J) m, n = size(A) lenI = length(I) @@ -2521,6 +2543,7 @@ _to_same_csc(::SparseMatrixCSC{Tv, Ti}, V::AbstractVector, I...) where {Tv,Ti} = setindex!(A::SparseMatrixCSC{Tv}, B::AbstractVecOrMat, I::Integer, J::Integer) where {Tv} = setindex!(A, convert(Tv, B), I, J) function setindex!(A::SparseMatrixCSC{Tv,Ti}, V::AbstractVecOrMat, Ix::Union{Integer, AbstractVector{<:Integer}, Colon}, Jx::Union{Integer, AbstractVector{<:Integer}, Colon}) where {Tv,Ti<:Integer} + @assert is_one_indexed(A, V, Ix, Jx) (I, J) = Base.ensure_indexable(to_indices(A, (Ix, Jx))) checkbounds(A, I, J) B = _to_same_csc(A, V, I, J) @@ -2654,6 +2677,7 @@ setindex!(A::Matrix, x::SparseMatrixCSC, I::AbstractVector{<:Integer}, J::Abstra setindex!(A::Matrix, x::SparseMatrixCSC, I::AbstractVector{Bool}, J::AbstractVector{<:Integer}) = setindex!(A, Array(x), findall(I), J) function setindex!(A::SparseMatrixCSC, x::AbstractArray, I::AbstractMatrix{Bool}) + @assert is_one_indexed(A, x, I) checkbounds(A, I) n = sum(I) (n == 0) && (return A) @@ -2754,6 +2778,7 @@ function setindex!(A::SparseMatrixCSC, x::AbstractArray, I::AbstractMatrix{Bool} end function setindex!(A::SparseMatrixCSC, x::AbstractArray, Ix::AbstractVector{<:Integer}) + @assert is_one_indexed(A, x, Ix) (I,) = Base.ensure_indexable(to_indices(A, (Ix,))) # We check bounds after sorting I n = length(I) @@ -2924,6 +2949,7 @@ julia> Base.SparseArrays.dropstored!(A, [1, 2], [1, 1]) """ function dropstored!(A::SparseMatrixCSC, I::AbstractVector{<:Integer}, J::AbstractVector{<:Integer}) + @assert is_one_indexed(A, I, J) m, n = size(A) nnzA = nnz(A) (nnzA == 0) && (return A) diff --git a/stdlib/SparseArrays/src/sparsevector.jl b/stdlib/SparseArrays/src/sparsevector.jl index f30ea679b88d2..327817c9d5c7b 100644 --- a/stdlib/SparseArrays/src/sparsevector.jl +++ b/stdlib/SparseArrays/src/sparsevector.jl @@ -182,6 +182,7 @@ julia> sparsevec([1, 3, 1, 2, 2], [true, true, false, false, false]) ``` """ function sparsevec(I::AbstractVector{<:Integer}, V::AbstractVector, combine::Function) + @assert is_one_indexed(I, V) length(I) == length(V) || throw(ArgumentError("index and value vectors must be the same length")) len = 0 @@ -195,6 +196,7 @@ function sparsevec(I::AbstractVector{<:Integer}, V::AbstractVector, combine::Fun end function sparsevec(I::AbstractVector{<:Integer}, V::AbstractVector, len::Integer, combine::Function) + @assert is_one_indexed(I, V) length(I) == length(V) || throw(ArgumentError("index and value vectors must be the same length")) for i in I @@ -377,6 +379,7 @@ sparsevec(a::AbstractSparseVector) = vec(a) sparse(a::AbstractVector) = sparsevec(a) function _dense2sparsevec(s::AbstractArray{Tv}, initcap::Ti) where {Tv,Ti} + @assert is_one_indexed(s) # pre-condition: initcap > 0; the initcap determines the index type n = length(s) cap = initcap @@ -534,6 +537,7 @@ end # Row slices getindex(A::SparseMatrixCSC, i::Integer, ::Colon) = A[i, 1:end] function Base.getindex(A::SparseMatrixCSC{Tv,Ti}, i::Integer, J::AbstractVector) where {Tv,Ti} + @assert is_one_indexed(A, J) checkbounds(A, i, J) nJ = length(J) colptrA = A.colptr; rowvalA = A.rowval; nzvalA = A.nzval @@ -567,6 +571,7 @@ end getindex(A::SparseMatrixCSC, I::AbstractVector{Bool}) = _logical_index(A, I) # Ambiguities getindex(A::SparseMatrixCSC, I::AbstractArray{Bool}) = _logical_index(A, I) function _logical_index(A::SparseMatrixCSC{Tv}, I::AbstractArray{Bool}) where Tv + @assert is_one_indexed(A, I) checkbounds(A, I) n = sum(I) nnzB = min(n, nnz(A)) @@ -608,6 +613,7 @@ end getindex(A::SparseMatrixCSC, ::Colon) = A[1:end] function getindex(A::SparseMatrixCSC{Tv}, I::AbstractUnitRange) where Tv + @assert is_one_indexed(A, I) checkbounds(A, I) szA = size(A) nA = szA[1]*szA[2] @@ -644,6 +650,7 @@ function getindex(A::SparseMatrixCSC{Tv}, I::AbstractUnitRange) where Tv end function getindex(A::SparseMatrixCSC{Tv,Ti}, I::AbstractVector) where {Tv,Ti} + @assert is_one_indexed(A, I) szA = size(A) nA = szA[1]*szA[2] colptrA = A.colptr @@ -857,6 +864,7 @@ end ### Conversion to matrix function SparseMatrixCSC{Tv,Ti}(x::AbstractSparseVector) where {Tv,Ti} + @assert is_one_indexed(x) n = length(x) xnzind = nonzeroinds(x) xnzval = nonzeros(x) @@ -874,6 +882,7 @@ SparseMatrixCSC{Tv}(x::AbstractSparseVector{<:Any,Ti}) where {Tv,Ti} = SparseMat SparseMatrixCSC(x::AbstractSparseVector{Tv,Ti}) where {Tv,Ti} = SparseMatrixCSC{Tv,Ti}(x) function Vector(x::AbstractSparseVector{Tv}) where Tv + @assert is_one_indexed(x) n = length(x) n == 0 && return Vector{Tv}() nzind = nonzeroinds(x) @@ -1075,6 +1084,7 @@ hvcat(rows::Tuple{Vararg{Int}}, xs::_TypedDenseConcatGroup{T}...) where {T} = Ba macro unarymap_nz2z_z2z(op, TF) esc(quote function $(op)(x::AbstractSparseVector{Tv,Ti}) where Tv<:$(TF) where Ti<:Integer + @assert is_one_indexed(x) R = typeof($(op)(zero(Tv))) xnzind = nonzeroinds(x) xnzval = nonzeros(x) @@ -1110,6 +1120,7 @@ imag(x::AbstractSparseVector{Tv,Ti}) where {Tv<:Real,Ti<:Integer} = SparseVector macro unarymap_z2nz(op, TF) esc(quote function $(op)(x::AbstractSparseVector{Tv,<:Integer}) where Tv<:$(TF) + @assert is_one_indexed(x) v0 = $(op)(zero(Tv)) R = typeof(v0) xnzind = nonzeroinds(x) @@ -1334,6 +1345,7 @@ adjoint(sv::SparseVector) = Adjoint(sv) # axpy function LinearAlgebra.axpy!(a::Number, x::SparseVectorUnion, y::AbstractVector) + @assert is_one_indexed(x, y) length(x) == length(y) || throw(DimensionMismatch()) nzind = nonzeroinds(x) nzval = nonzeros(x) @@ -1387,6 +1399,7 @@ end # dot function dot(x::AbstractVector{Tx}, y::SparseVectorUnion{Ty}) where {Tx<:Number,Ty<:Number} + @assert is_one_indexed(x, y) n = length(x) length(y) == n || throw(DimensionMismatch()) nzind = nonzeroinds(y) @@ -1399,6 +1412,7 @@ function dot(x::AbstractVector{Tx}, y::SparseVectorUnion{Ty}) where {Tx<:Number, end function dot(x::SparseVectorUnion{Tx}, y::AbstractVector{Ty}) where {Tx<:Number,Ty<:Number} + @assert is_one_indexed(x, y) n = length(y) length(x) == n || throw(DimensionMismatch()) nzind = nonzeroinds(x) @@ -1451,6 +1465,7 @@ end # lowrankupdate (BLAS.ger! like) function LinearAlgebra.lowrankupdate!(A::StridedMatrix, x::AbstractVector, y::SparseVectorUnion, α::Number = 1) + @assert is_one_indexed(A, x, y) nzi = nonzeroinds(y) nzv = nonzeros(y) @inbounds for (j,v) in zip(nzi,nzv) @@ -1465,6 +1480,7 @@ end # * and mul! function (*)(A::StridedMatrix{Ta}, x::AbstractSparseVector{Tx}) where {Ta,Tx} + @assert is_one_indexed(A, x) m, n = size(A) length(x) == n || throw(DimensionMismatch()) Ty = promote_type(Ta, Tx) @@ -1476,6 +1492,7 @@ mul!(y::AbstractVector{Ty}, A::StridedMatrix, x::AbstractSparseVector{Tx}) where mul!(y, A, x, one(Tx), zero(Ty)) function mul!(y::AbstractVector, A::StridedMatrix, x::AbstractSparseVector, α::Number, β::Number) + @assert is_one_indexed(y, A, x) m, n = size(A) length(x) == n && length(y) == m || throw(DimensionMismatch()) m == 0 && return y @@ -1502,6 +1519,7 @@ end # * and mul!(C, transpose(A), B) function *(transA::Transpose{<:Any,<:StridedMatrix{Ta}}, x::AbstractSparseVector{Tx}) where {Ta,Tx} + @assert is_one_indexed(transA, x) A = transA.parent m, n = size(A) length(x) == m || throw(DimensionMismatch()) @@ -1515,6 +1533,7 @@ mul!(y::AbstractVector{Ty}, transA::Transpose{<:Any,<:StridedMatrix}, x::Abstrac function mul!(y::AbstractVector, transA::Transpose{<:Any,<:StridedMatrix}, x::AbstractSparseVector, α::Number, β::Number) A = transA.parent + @assert is_one_indexed(y, A, x) m, n = size(A) length(x) == m && length(y) == n || throw(DimensionMismatch()) n == 0 && return y @@ -1544,6 +1563,7 @@ end function densemv(A::SparseMatrixCSC, x::AbstractSparseVector; trans::AbstractChar='N') local xlen::Int, ylen::Int + @assert is_one_indexed(A, x) m, n = size(A) if trans == 'N' || trans == 'n' xlen = n; ylen = m @@ -1573,6 +1593,7 @@ mul!(y::AbstractVector{Ty}, A::SparseMatrixCSC, x::AbstractSparseVector{Tx}) whe mul!(y, A, x, one(Tx), zero(Ty)) function mul!(y::AbstractVector, A::SparseMatrixCSC, x::AbstractSparseVector, α::Number, β::Number) + @assert is_one_indexed(y, A, x) m, n = size(A) length(x) == n && length(y) == m || throw(DimensionMismatch()) m == 0 && return y @@ -1617,6 +1638,7 @@ mul!(y::AbstractVector, adjA::Adjoint{<:Any,<:SparseMatrixCSC}, x::AbstractSpars function _At_or_Ac_mul_B!(tfun::Function, y::AbstractVector, A::SparseMatrixCSC, x::AbstractSparseVector, α::Number, β::Number) + @assert is_one_indexed(y, A, x) m, n = size(A) length(x) == m && length(y) == n || throw(DimensionMismatch()) n == 0 && return y @@ -1645,6 +1667,7 @@ end ### BLAS-2 / sparse A * sparse x -> dense y function *(A::SparseMatrixCSC, x::AbstractSparseVector) + @assert is_one_indexed(A, x) y = densemv(A, x) initcap = min(nnz(A), size(A,1)) _dense2sparsevec(y, initcap) @@ -1657,6 +1680,7 @@ end (A = adjA.parent; _At_or_Ac_mul_B(dot, A, x)) function _At_or_Ac_mul_B(tfun::Function, A::SparseMatrixCSC{TvA,TiA}, x::AbstractSparseVector{TvX,TiX}) where {TvA,TiA,TvX,TiX} + @assert is_one_indexed(A, x) m, n = size(A) length(x) == m || throw(DimensionMismatch()) Tv = promote_type(TvA, TvX) diff --git a/stdlib/Statistics/src/Statistics.jl b/stdlib/Statistics/src/Statistics.jl index 85d8cd2bcd9f0..c9b43620a3572 100644 --- a/stdlib/Statistics/src/Statistics.jl +++ b/stdlib/Statistics/src/Statistics.jl @@ -196,7 +196,7 @@ function centralize_sumabs2!(R::AbstractArray{S}, A::AbstractArray, means::Abstr isempty(R) || fill!(R, zero(S)) isempty(A) && return R - if Base.has_fast_linear_indexing(A) && lsiz > 16 + if Base.has_fast_linear_indexing(A) && lsiz > 16 && is_one_indexed(R, means) nslices = div(Base._length(A), lsiz) ibase = first(LinearIndices(A))-1 for i = 1:nslices @@ -506,6 +506,7 @@ clampcor(x) = x # cov2cor! function cov2cor!(C::AbstractMatrix{T}, xsd::AbstractArray) where T + @assert is_one_indexed(C, xsd) nx = length(xsd) size(C) == (nx, nx) || throw(DimensionMismatch("inconsistent dimensions")) for j = 1:nx @@ -520,6 +521,7 @@ function cov2cor!(C::AbstractMatrix{T}, xsd::AbstractArray) where T return C end function cov2cor!(C::AbstractMatrix, xsd, ysd::AbstractArray) + @assert is_one_indexed(C, ysd) nx, ny = size(C) length(ysd) == ny || throw(DimensionMismatch("inconsistent dimensions")) for (j, y) in enumerate(ysd) # fixme (iter): here and in all `cov2cor!` we assume that `C` is efficiently indexed by integers @@ -530,6 +532,7 @@ function cov2cor!(C::AbstractMatrix, xsd, ysd::AbstractArray) return C end function cov2cor!(C::AbstractMatrix, xsd::AbstractArray, ysd) + @assert is_one_indexed(C, xsd) nx, ny = size(C) length(xsd) == nx || throw(DimensionMismatch("inconsistent dimensions")) for j in 1:ny @@ -540,6 +543,7 @@ function cov2cor!(C::AbstractMatrix, xsd::AbstractArray, ysd) return C end function cov2cor!(C::AbstractMatrix, xsd::AbstractArray, ysd::AbstractArray) + @assert is_one_indexed(C, xsd, ysd) nx, ny = size(C) (length(xsd) == nx && length(ysd) == ny) || throw(DimensionMismatch("inconsistent dimensions")) @@ -570,6 +574,7 @@ corzm(x::AbstractMatrix, y::AbstractMatrix, vardim::Int=1) = corm(x::AbstractVector{T}, xmean) where {T} = one(real(T)) corm(x::AbstractMatrix, xmean, vardim::Int=1) = corzm(x .- xmean, vardim) function corm(x::AbstractVector, mx, y::AbstractVector, my) + @assert is_one_indexed(x, y) n = length(x) length(y) == n || throw(DimensionMismatch("inconsistent lengths")) n > 0 || throw(ArgumentError("correlation only defined for non-empty vectors")) @@ -794,6 +799,7 @@ julia> y """ function quantile!(q::AbstractArray, v::AbstractVector, p::AbstractArray; sorted::Bool=false) + @assert is_one_indexed(q, v, p) if size(p) != size(q) throw(DimensionMismatch("size of p, $(size(p)), must equal size of q, $(size(q))")) end @@ -824,6 +830,7 @@ end # Function to perform partial sort of v for quantiles in given range function _quantilesort!(v::AbstractArray, sorted::Bool, minp::Real, maxp::Real) isempty(v) && throw(ArgumentError("empty data vector")) + @assert is_one_indexed(v) if !sorted lv = length(v) @@ -841,6 +848,7 @@ end # Core quantile lookup function: assumes `v` sorted @inline function _quantile(v::AbstractVector, p::Real) 0 <= p <= 1 || throw(ArgumentError("input probability out of [0,1] range")) + @assert is_one_indexed(v) lv = length(v) f0 = (lv - 1)*p # 0-based interpolated index @@ -932,6 +940,7 @@ end # This is the function that does the reduction underlying var/std function centralize_sumabs2!(R::AbstractArray{S}, A::SparseMatrixCSC{Tv,Ti}, means::AbstractArray) where {S,Tv,Ti} + @assert is_one_indexed(R, A, means) lsiz = Base.check_reducedims(R,A) size(means) == size(R) || error("size of means must match size of R") isempty(R) || fill!(R, zero(S)) diff --git a/stdlib/SuiteSparse/src/SuiteSparse.jl b/stdlib/SuiteSparse/src/SuiteSparse.jl index 99a0128d73f47..77e2a8461bfae 100644 --- a/stdlib/SuiteSparse/src/SuiteSparse.jl +++ b/stdlib/SuiteSparse/src/SuiteSparse.jl @@ -11,14 +11,14 @@ import LinearAlgebra: ldiv!, rdiv! # Convert from 1-based to 0-based indices function decrement!(A::AbstractArray{T}) where T<:Integer - for i in 1:length(A); A[i] -= oneunit(T) end + for i in eachindex(A); A[i] -= oneunit(T) end A end decrement(A::AbstractArray{<:Integer}) = decrement!(copy(A)) # Convert from 0-based to 1-based indices function increment!(A::AbstractArray{T}) where T<:Integer - for i in 1:length(A); A[i] += oneunit(T) end + for i in eachindex(A); A[i] += oneunit(T) end A end increment(A::AbstractArray{<:Integer}) = increment!(copy(A)) diff --git a/stdlib/SuiteSparse/src/cholmod.jl b/stdlib/SuiteSparse/src/cholmod.jl index c64f199cdc818..3edce415c0183 100644 --- a/stdlib/SuiteSparse/src/cholmod.jl +++ b/stdlib/SuiteSparse/src/cholmod.jl @@ -1040,6 +1040,7 @@ Base.copyto!(dest::AbstractArray{T,2}, D::Dense{T}) where {T<:VTypes} = _copy!(d Base.copyto!(dest::AbstractArray, D::Dense) = _copy!(dest, D) function _copy!(dest::AbstractArray, D::Dense) + @assert is_one_indexed(dest) s = unsafe_load(pointer(D)) n = s.nrow*s.ncol n <= length(dest) || throw(BoundsError(dest, n)) diff --git a/stdlib/SuiteSparse/src/spqr.jl b/stdlib/SuiteSparse/src/spqr.jl index d2a7b923591f1..5b5665b3386e8 100644 --- a/stdlib/SuiteSparse/src/spqr.jl +++ b/stdlib/SuiteSparse/src/spqr.jl @@ -346,6 +346,7 @@ function (\)(F::QRSparse{T}, B::VecOrMat{Complex{T}}) where T<:LinearAlgebra.Bla # |z2|z4| -> |y1|y2|y3|y4| -> |x2|y2| -> |x2|y2|x4|y4| # |x3|y3| # |x4|y4| + @assert is_one_indexed(F, B) c2r = reshape(copy(transpose(reinterpret(T, reshape(B, (1, length(B)))))), size(B, 1), 2*size(B, 2)) x = F\c2r