From c168e712b0b6eef013a1e315fc69cf863f479141 Mon Sep 17 00:00:00 2001 From: Andy Ferris Date: Wed, 21 Dec 2016 16:11:48 +1000 Subject: [PATCH 1/2] Moved all array transpose functions to LinAlg Transposition is a concept of linear algebra rather than multidimensional arrays of data. --- base/abstractarray.jl | 23 ------ base/abstractarraymath.jl | 2 - base/arraymath.jl | 127 ------------------------------- base/bitarray.jl | 88 --------------------- base/linalg/bitarray.jl | 89 ++++++++++++++++++++++ base/linalg/linalg.jl | 15 ++-- base/linalg/transpose.jl | 156 ++++++++++++++++++++++++++++++++++++++ base/range.jl | 3 - 8 files changed, 255 insertions(+), 248 deletions(-) create mode 100644 base/linalg/transpose.jl diff --git a/base/abstractarray.jl b/base/abstractarray.jl index 64d58c2bc8a95..f678b471e817e 100644 --- a/base/abstractarray.jl +++ b/base/abstractarray.jl @@ -643,29 +643,6 @@ function copy!{R,S}(B::AbstractVecOrMat{R}, ir_dest::Range{Int}, jr_dest::Range{ return B end -function copy_transpose!{R,S}(B::AbstractVecOrMat{R}, ir_dest::Range{Int}, jr_dest::Range{Int}, - A::AbstractVecOrMat{S}, ir_src::Range{Int}, jr_src::Range{Int}) - if length(ir_dest) != length(jr_src) - throw(ArgumentError(string("source and destination must have same size (got ", - length(jr_src)," and ",length(ir_dest),")"))) - end - if length(jr_dest) != length(ir_src) - throw(ArgumentError(string("source and destination must have same size (got ", - length(ir_src)," and ",length(jr_dest),")"))) - end - @boundscheck checkbounds(B, ir_dest, jr_dest) - @boundscheck checkbounds(A, ir_src, jr_src) - idest = first(ir_dest) - for jsrc in jr_src - jdest = first(jr_dest) - for isrc in ir_src - B[idest,jdest] = A[isrc,jsrc] - jdest += step(jr_dest) - end - idest += step(ir_dest) - end - return B -end """ copymutable(a) diff --git a/base/abstractarraymath.jl b/base/abstractarraymath.jl index 6580333dd0a65..93135aed45faa 100644 --- a/base/abstractarraymath.jl +++ b/base/abstractarraymath.jl @@ -7,8 +7,6 @@ isinteger{T<:Integer,n}(x::AbstractArray{T,n}) = true isreal(x::AbstractArray) = all(isreal,x) iszero(x::AbstractArray) = all(iszero,x) isreal{T<:Real,n}(x::AbstractArray{T,n}) = true -ctranspose(a::AbstractArray) = error("ctranspose not implemented for $(typeof(a)). Consider adding parentheses, e.g. A*(B*C') instead of A*B*C' to avoid explicit calculation of the transposed matrix.") -transpose(a::AbstractArray) = error("transpose not implemented for $(typeof(a)). Consider adding parentheses, e.g. A*(B*C.') instead of A*B*C' to avoid explicit calculation of the transposed matrix.") ## Constructors ## diff --git a/base/arraymath.jl b/base/arraymath.jl index 82685dd6a69b1..1e47fff0ddbb1 100644 --- a/base/arraymath.jl +++ b/base/arraymath.jl @@ -274,130 +274,3 @@ julia> rot180(a,2) ``` """ rot180(A::AbstractMatrix, k::Integer) = mod(k, 2) == 1 ? rot180(A) : copy(A) - -## Transpose ## - -""" - transpose!(dest,src) - -Transpose array `src` and store the result in the preallocated array `dest`, which should -have a size corresponding to `(size(src,2),size(src,1))`. No in-place transposition is -supported and unexpected results will happen if `src` and `dest` have overlapping memory -regions. -""" -transpose!(B::AbstractMatrix, A::AbstractMatrix) = transpose_f!(transpose, B, A) - -""" - ctranspose!(dest,src) - -Conjugate transpose array `src` and store the result in the preallocated array `dest`, which -should have a size corresponding to `(size(src,2),size(src,1))`. No in-place transposition -is supported and unexpected results will happen if `src` and `dest` have overlapping memory -regions. -""" -ctranspose!(B::AbstractMatrix, A::AbstractMatrix) = transpose_f!(ctranspose, B, A) -function transpose!(B::AbstractVector, A::AbstractMatrix) - indices(B,1) == indices(A,2) && indices(A,1) == 1:1 || throw(DimensionMismatch("transpose")) - copy!(B, A) -end -function transpose!(B::AbstractMatrix, A::AbstractVector) - indices(B,2) == indices(A,1) && indices(B,1) == 1:1 || throw(DimensionMismatch("transpose")) - copy!(B, A) -end -function ctranspose!(B::AbstractVector, A::AbstractMatrix) - indices(B,1) == indices(A,2) && indices(A,1) == 1:1 || throw(DimensionMismatch("transpose")) - ccopy!(B, A) -end -function ctranspose!(B::AbstractMatrix, A::AbstractVector) - indices(B,2) == indices(A,1) && indices(B,1) == 1:1 || throw(DimensionMismatch("transpose")) - ccopy!(B, A) -end - -const transposebaselength=64 -function transpose_f!(f,B::AbstractMatrix,A::AbstractMatrix) - inds = indices(A) - indices(B,1) == inds[2] && indices(B,2) == inds[1] || throw(DimensionMismatch(string(f))) - - m, n = length(inds[1]), length(inds[2]) - if m*n<=4*transposebaselength - @inbounds begin - for j = inds[2] - for i = inds[1] - B[j,i] = f(A[i,j]) - end - end - end - else - transposeblock!(f,B,A,m,n,first(inds[1])-1,first(inds[2])-1) - end - return B -end -function transposeblock!(f,B::AbstractMatrix,A::AbstractMatrix,m::Int,n::Int,offseti::Int,offsetj::Int) - if m*n<=transposebaselength - @inbounds begin - for j = offsetj+(1:n) - for i = offseti+(1:m) - B[j,i] = f(A[i,j]) - end - end - end - elseif m>n - newm=m>>1 - transposeblock!(f,B,A,newm,n,offseti,offsetj) - transposeblock!(f,B,A,m-newm,n,offseti+newm,offsetj) - else - newn=n>>1 - transposeblock!(f,B,A,m,newn,offseti,offsetj) - transposeblock!(f,B,A,m,n-newn,offseti,offsetj+newn) - end - return B -end - -function ccopy!(B, A) - RB, RA = eachindex(B), eachindex(A) - if RB == RA - for i = RB - B[i] = ctranspose(A[i]) - end - else - for (i,j) = zip(RB, RA) - B[i] = ctranspose(A[j]) - end - end -end - -""" - transpose(A) - -The transposition operator (`.'`). - -# Example - -```jldoctest -julia> A = [1 2 3; 4 5 6; 7 8 9] -3×3 Array{Int64,2}: - 1 2 3 - 4 5 6 - 7 8 9 - -julia> transpose(A) -3×3 Array{Int64,2}: - 1 4 7 - 2 5 8 - 3 6 9 -``` -""" -function transpose(A::AbstractMatrix) - ind1, ind2 = indices(A) - B = similar(A, (ind2, ind1)) - transpose!(B, A) -end -function ctranspose(A::AbstractMatrix) - ind1, ind2 = indices(A) - B = similar(A, (ind2, ind1)) - ctranspose!(B, A) -end -ctranspose{T<:Real}(A::AbstractVecOrMat{T}) = transpose(A) - -transpose(x::AbstractVector) = [ transpose(v) for i=of_indices(x, OneTo(1)), v in x ] -ctranspose{T}(x::AbstractVector{T}) = T[ ctranspose(v) for i=of_indices(x, OneTo(1)), v in x ] diff --git a/base/bitarray.jl b/base/bitarray.jl index d41fa8f64308a..1a0c6b4940a89 100644 --- a/base/bitarray.jl +++ b/base/bitarray.jl @@ -1893,94 +1893,6 @@ function filter(f, Bs::BitArray) Bs[boolmap] end -## Transpose ## - -transpose(B::BitVector) = reshape(copy(B), 1, length(B)) - -# fast 8x8 bit transpose from Henry S. Warrens's "Hacker's Delight" -# http://www.hackersdelight.org/hdcodetxt/transpose8.c.txt -function transpose8x8(x::UInt64) - y = x - t = xor(y, y >>> 7) & 0x00aa00aa00aa00aa - y = xor(y, t, t << 7) - t = xor(y, y >>> 14) & 0x0000cccc0000cccc - y = xor(y, t, t << 14) - t = xor(y, y >>> 28) & 0x00000000f0f0f0f0 - return xor(y, t, t << 28) -end - -function form_8x8_chunk(Bc::Vector{UInt64}, i1::Int, i2::Int, m::Int, cgap::Int, cinc::Int, nc::Int, msk8::UInt64) - x = UInt64(0) - - k, l = get_chunks_id(i1 + (i2 - 1) * m) - r = 0 - for j = 1:8 - k > nc && break - x |= ((Bc[k] >>> l) & msk8) << r - if l + 8 >= 64 && nc > k - r0 = 8 - _mod64(l + 8) - x |= (Bc[k + 1] & (msk8 >>> r0)) << (r + r0) - end - k += cgap + (l + cinc >= 64 ? 1 : 0) - l = _mod64(l + cinc) - r += 8 - end - return x -end - -# note: assumes B is filled with 0's -function put_8x8_chunk(Bc::Vector{UInt64}, i1::Int, i2::Int, x::UInt64, m::Int, cgap::Int, cinc::Int, nc::Int, msk8::UInt64) - k, l = get_chunks_id(i1 + (i2 - 1) * m) - r = 0 - for j = 1:8 - k > nc && break - Bc[k] |= ((x >>> r) & msk8) << l - if l + 8 >= 64 && nc > k - r0 = 8 - _mod64(l + 8) - Bc[k + 1] |= ((x >>> (r + r0)) & (msk8 >>> r0)) - end - k += cgap + (l + cinc >= 64 ? 1 : 0) - l = _mod64(l + cinc) - r += 8 - end - return -end - -function transpose(B::BitMatrix) - l1 = size(B, 1) - l2 = size(B, 2) - Bt = falses(l2, l1) - - cgap1, cinc1 = _div64(l1), _mod64(l1) - cgap2, cinc2 = _div64(l2), _mod64(l2) - - Bc = B.chunks - Btc = Bt.chunks - - nc = length(Bc) - - for i = 1:8:l1 - msk8_1 = UInt64(0xff) - if (l1 < i + 7) - msk8_1 >>>= i + 7 - l1 - end - - for j = 1:8:l2 - x = form_8x8_chunk(Bc, i, j, l1, cgap1, cinc1, nc, msk8_1) - x = transpose8x8(x) - - msk8_2 = UInt64(0xff) - if (l2 < j + 7) - msk8_2 >>>= j + 7 - l2 - end - - put_8x8_chunk(Btc, j, i, x, l2, cgap2, cinc2, nc, msk8_2) - end - end - return Bt -end - -ctranspose(B::BitArray) = transpose(B) ## Concatenation ## diff --git a/base/linalg/bitarray.jl b/base/linalg/bitarray.jl index cb59cf4a830a4..72124add2aac0 100644 --- a/base/linalg/bitarray.jl +++ b/base/linalg/bitarray.jl @@ -211,3 +211,92 @@ function findmin(a::BitArray) k == l || return (false, ti) return m, mi end + +## transpose ## + +transpose(B::BitVector) = reshape(copy(B), 1, length(B)) + +# fast 8x8 bit transpose from Henry S. Warrens's "Hacker's Delight" +# http://www.hackersdelight.org/hdcodetxt/transpose8.c.txt +function transpose8x8(x::UInt64) + y = x + t = xor(y, y >>> 7) & 0x00aa00aa00aa00aa + y = xor(y, t, t << 7) + t = xor(y, y >>> 14) & 0x0000cccc0000cccc + y = xor(y, t, t << 14) + t = xor(y, y >>> 28) & 0x00000000f0f0f0f0 + return xor(y, t, t << 28) +end + +function form_8x8_chunk(Bc::Vector{UInt64}, i1::Int, i2::Int, m::Int, cgap::Int, cinc::Int, nc::Int, msk8::UInt64) + x = UInt64(0) + + k, l = Base.get_chunks_id(i1 + (i2 - 1) * m) + r = 0 + for j = 1:8 + k > nc && break + x |= ((Bc[k] >>> l) & msk8) << r + if l + 8 >= 64 && nc > k + r0 = 8 - Base._mod64(l + 8) + x |= (Bc[k + 1] & (msk8 >>> r0)) << (r + r0) + end + k += cgap + (l + cinc >= 64 ? 1 : 0) + l = Base._mod64(l + cinc) + r += 8 + end + return x +end + +# note: assumes B is filled with 0's +function put_8x8_chunk(Bc::Vector{UInt64}, i1::Int, i2::Int, x::UInt64, m::Int, cgap::Int, cinc::Int, nc::Int, msk8::UInt64) + k, l = Base.get_chunks_id(i1 + (i2 - 1) * m) + r = 0 + for j = 1:8 + k > nc && break + Bc[k] |= ((x >>> r) & msk8) << l + if l + 8 >= 64 && nc > k + r0 = 8 - Base._mod64(l + 8) + Bc[k + 1] |= ((x >>> (r + r0)) & (msk8 >>> r0)) + end + k += cgap + (l + cinc >= 64 ? 1 : 0) + l = Base._mod64(l + cinc) + r += 8 + end + return +end + +function transpose(B::BitMatrix) + l1 = size(B, 1) + l2 = size(B, 2) + Bt = falses(l2, l1) + + cgap1, cinc1 = Base._div64(l1), Base._mod64(l1) + cgap2, cinc2 = Base._div64(l2), Base._mod64(l2) + + Bc = B.chunks + Btc = Bt.chunks + + nc = length(Bc) + + for i = 1:8:l1 + msk8_1 = UInt64(0xff) + if (l1 < i + 7) + msk8_1 >>>= i + 7 - l1 + end + + for j = 1:8:l2 + x = form_8x8_chunk(Bc, i, j, l1, cgap1, cinc1, nc, msk8_1) + x = transpose8x8(x) + + msk8_2 = UInt64(0xff) + if (l2 < j + 7) + msk8_2 >>>= j + 7 - l2 + end + + put_8x8_chunk(Btc, j, i, x, l2, cgap2, cinc2, nc, msk8_2) + end + end + return Bt +end + +ctranspose(B::BitArray) = transpose(B) diff --git a/base/linalg/linalg.jl b/base/linalg/linalg.jl index 2fc634e64e559..e817b3a83d118 100644 --- a/base/linalg/linalg.jl +++ b/base/linalg/linalg.jl @@ -5,11 +5,11 @@ module LinAlg import Base: \, /, *, ^, +, -, == import Base: A_mul_Bt, At_ldiv_Bt, A_rdiv_Bc, At_ldiv_B, Ac_mul_Bc, A_mul_Bc, Ac_mul_B, Ac_ldiv_B, Ac_ldiv_Bc, At_mul_Bt, A_rdiv_Bt, At_mul_B -import Base: USE_BLAS64, abs, big, ceil, conj, convert, copy, copy!, copy_transpose!, - ctranspose, ctranspose!, eltype, eye, findmax, findmin, fill!, floor, full, getindex, - imag, inv, isapprox, kron, ndims, parent, power_by_squaring, print_matrix, - promote_rule, real, round, setindex!, show, similar, size, transpose, transpose!, - trunc, broadcast +import Base: USE_BLAS64, abs, big, ceil, conj, convert, copy, copy!, ctranspose, + eltype, eye, findmax, findmin, fill!, floor, full, getindex, imag, inv, + isapprox, kron, ndims, parent, power_by_squaring, print_matrix, + promote_rule, real, round, setindex!, show, similar, size, transpose, trunc, + broadcast using Base: promote_op, _length, iszero # 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))`. @@ -56,8 +56,10 @@ export cond, condskeel, copy!, + copy_transpose!, cross, ctranspose, + ctranspose!, det, diag, diagind, @@ -127,6 +129,7 @@ export sylvester, trace, transpose, + transpose!, tril, triu, tril!, @@ -231,6 +234,8 @@ end copy_oftype{T,N}(A::AbstractArray{T,N}, ::Type{T}) = copy(A) copy_oftype{T,N,S}(A::AbstractArray{T,N}, ::Type{S}) = convert(AbstractArray{S,N}, A) +include("transpose.jl") + include("exceptions.jl") include("generic.jl") diff --git a/base/linalg/transpose.jl b/base/linalg/transpose.jl new file mode 100644 index 0000000000000..ab4007486509b --- /dev/null +++ b/base/linalg/transpose.jl @@ -0,0 +1,156 @@ +# This file is a part of Julia. License is MIT: http://julialang.org/license + +ctranspose(a::AbstractArray) = error("ctranspose not defined for $(typeof(a)). Consider using `permutedims` for higher-dimensional arrays.") +transpose(a::AbstractArray) = error("transpose not defined for $(typeof(a)). Consider using `permutedims` for higher-dimensional arrays.") + +## Matrix transposition ## + +""" + transpose!(dest,src) + +Transpose array `src` and store the result in the preallocated array `dest`, which should +have a size corresponding to `(size(src,2),size(src,1))`. No in-place transposition is +supported and unexpected results will happen if `src` and `dest` have overlapping memory +regions. +""" +transpose!(B::AbstractMatrix, A::AbstractMatrix) = transpose_f!(transpose, B, A) + +""" + ctranspose!(dest,src) + +Conjugate transpose array `src` and store the result in the preallocated array `dest`, which +should have a size corresponding to `(size(src,2),size(src,1))`. No in-place transposition +is supported and unexpected results will happen if `src` and `dest` have overlapping memory +regions. +""" +ctranspose!(B::AbstractMatrix, A::AbstractMatrix) = transpose_f!(ctranspose, B, A) +function transpose!(B::AbstractVector, A::AbstractMatrix) + indices(B,1) == indices(A,2) && indices(A,1) == 1:1 || throw(DimensionMismatch("transpose")) + copy!(B, A) +end +function transpose!(B::AbstractMatrix, A::AbstractVector) + indices(B,2) == indices(A,1) && indices(B,1) == 1:1 || throw(DimensionMismatch("transpose")) + copy!(B, A) +end +function ctranspose!(B::AbstractVector, A::AbstractMatrix) + indices(B,1) == indices(A,2) && indices(A,1) == 1:1 || throw(DimensionMismatch("transpose")) + ccopy!(B, A) +end +function ctranspose!(B::AbstractMatrix, A::AbstractVector) + indices(B,2) == indices(A,1) && indices(B,1) == 1:1 || throw(DimensionMismatch("transpose")) + ccopy!(B, A) +end + +const transposebaselength=64 +function transpose_f!(f,B::AbstractMatrix,A::AbstractMatrix) + inds = indices(A) + indices(B,1) == inds[2] && indices(B,2) == inds[1] || throw(DimensionMismatch(string(f))) + + m, n = length(inds[1]), length(inds[2]) + if m*n<=4*transposebaselength + @inbounds begin + for j = inds[2] + for i = inds[1] + B[j,i] = f(A[i,j]) + end + end + end + else + transposeblock!(f,B,A,m,n,first(inds[1])-1,first(inds[2])-1) + end + return B +end +function transposeblock!(f,B::AbstractMatrix,A::AbstractMatrix,m::Int,n::Int,offseti::Int,offsetj::Int) + if m*n<=transposebaselength + @inbounds begin + for j = offsetj+(1:n) + for i = offseti+(1:m) + B[j,i] = f(A[i,j]) + end + end + end + elseif m>n + newm=m>>1 + transposeblock!(f,B,A,newm,n,offseti,offsetj) + transposeblock!(f,B,A,m-newm,n,offseti+newm,offsetj) + else + newn=n>>1 + transposeblock!(f,B,A,m,newn,offseti,offsetj) + transposeblock!(f,B,A,m,n-newn,offseti,offsetj+newn) + end + return B +end + +function ccopy!(B, A) + RB, RA = eachindex(B), eachindex(A) + if RB == RA + for i = RB + B[i] = ctranspose(A[i]) + end + else + for (i,j) = zip(RB, RA) + B[i] = ctranspose(A[j]) + end + end +end + +""" + transpose(A::AbstractMatrix) + +The transposition operator (`.'`). + +# Example + +```jldoctest +julia> A = [1 2 3; 4 5 6; 7 8 9] +3×3 Array{Int64,2}: + 1 2 3 + 4 5 6 + 7 8 9 + +julia> transpose(A) +3×3 Array{Int64,2}: + 1 4 7 + 2 5 8 + 3 6 9 +``` +""" +function transpose(A::AbstractMatrix) + ind1, ind2 = indices(A) + B = similar(A, (ind2, ind1)) + transpose!(B, A) +end +function ctranspose(A::AbstractMatrix) + ind1, ind2 = indices(A) + B = similar(A, (ind2, ind1)) + ctranspose!(B, A) +end + +@inline ctranspose{T<:Real}(A::AbstractVecOrMat{T}) = transpose(A) + +transpose(x::AbstractVector) = [ transpose(v) for i=Base.of_indices(x, Base.OneTo(1)), v in x ] +ctranspose{T}(x::AbstractVector{T}) = T[ ctranspose(v) for i=Base.of_indices(x, Base.OneTo(1)), v in x ] + +function copy_transpose!{R,S}(B::AbstractVecOrMat{R}, ir_dest::Range{Int}, jr_dest::Range{Int}, + A::AbstractVecOrMat{S}, ir_src::Range{Int}, jr_src::Range{Int}) + if length(ir_dest) != length(jr_src) + throw(ArgumentError(string("source and destination must have same size (got ", + length(jr_src)," and ",length(ir_dest),")"))) + end + if length(jr_dest) != length(ir_src) + throw(ArgumentError(string("source and destination must have same size (got ", + length(ir_src)," and ",length(jr_dest),")"))) + end + @boundscheck checkbounds(B, ir_dest, jr_dest) + @boundscheck checkbounds(A, ir_src, jr_src) + idest = first(ir_dest) + for jsrc in jr_src + jdest = first(jr_dest) + for isrc in ir_src + B[idest,jdest] = A[isrc,jsrc] + jdest += step(jr_dest) + end + idest += step(ir_dest) + end + return B +end diff --git a/base/range.jl b/base/range.jl index 48aaf0de11592..d58238d93ee0f 100644 --- a/base/range.jl +++ b/base/range.jl @@ -463,9 +463,6 @@ maximum(r::AbstractUnitRange) = isempty(r) ? throw(ArgumentError("range must be minimum(r::Range) = isempty(r) ? throw(ArgumentError("range must be non-empty")) : min(first(r), last(r)) maximum(r::Range) = isempty(r) ? throw(ArgumentError("range must be non-empty")) : max(first(r), last(r)) -ctranspose(r::Range) = [x for _=1:1, x=r] -transpose(r::Range) = r' - # Ranges are immutable copy(r::Range) = r From af9a28f761f78ea1226f0f2bd9115f53732919c8 Mon Sep 17 00:00:00 2001 From: Andy Ferris Date: Wed, 21 Dec 2016 17:08:32 +1000 Subject: [PATCH 2/2] Introduce `RowVector` as vector transpose `RowVector` is now defined as the `transpose` of any `AbstractVector`. If `v` is an `AbstractVector`, then it obeys the identity that `(v.').' === v` and the matrix multiplication rules follow that `(A * v).' == (v.' * A.')`. `RowVector` is a "view" and maintains the recursive nature of `transpose`. It is a subtype of `AbstractMatrix` and maintains the current shape and broadcast behavior for `v.'. Consequences include: * v'' is a vector, not a matrix * v'*v is a scalar, not a vector * v*v' is the outer produce (returns a matrix) * v*A (for A::AbstractMatrix) is removed, since its puprose was to provide a way of doing the outer product above, and is no longer necessary. Closes #4774 --- NEWS.md | 4 + base/docs/basedocs.jl | 10 ++ base/exports.jl | 1 + base/linalg/bidiag.jl | 9 ++ base/linalg/bitarray.jl | 4 - base/linalg/dense.jl | 2 +- base/linalg/diagonal.jl | 12 ++ base/linalg/generic.jl | 14 +++ base/linalg/linalg.jl | 16 ++- base/linalg/matmul.jl | 5 +- base/linalg/qr.jl | 2 + base/linalg/rowvector.jl | 214 +++++++++++++++++++++++++++++++ base/linalg/transpose.jl | 3 - base/linalg/triangular.jl | 41 ++++++ base/sparse/cholmod.jl | 1 + base/sparse/linalg.jl | 2 + base/sparse/sparsevector.jl | 14 +-- doc/src/stdlib/linalg.md | 1 + test/arrayops.jl | 6 +- test/bitarray.jl | 5 +- test/choosetests.jl | 2 +- test/linalg/dense.jl | 2 +- test/linalg/matmul.jl | 2 +- test/linalg/rowvector.jl | 245 ++++++++++++++++++++++++++++++++++++ test/offsetarray.jl | 9 +- test/show.jl | 2 +- test/sparse/sparse.jl | 9 +- 27 files changed, 590 insertions(+), 47 deletions(-) create mode 100644 base/linalg/rowvector.jl create mode 100644 test/linalg/rowvector.jl diff --git a/NEWS.md b/NEWS.md index 48c35498e834a..6580b32e3017d 100644 --- a/NEWS.md +++ b/NEWS.md @@ -99,6 +99,10 @@ This section lists changes that do not have deprecation warnings. `n` and `max_delay`. The previous functionality can be achieved setting `delays` to `ExponentialBackOff`. ([#19331]) + * `transpose(::AbstractVector)` now always returns a `RowVector` view of the input (which is a + special 1×n-sized `AbstractMatrix`), not a `Matrix`, etc. In particular, for + `v::AbstractVector` we now have `(v.').' === v` and `v.' * v` is a scalar. ([#19670]) + Library improvements -------------------- diff --git a/base/docs/basedocs.jl b/base/docs/basedocs.jl index f1198a74956cf..cc75293a2ca02 100644 --- a/base/docs/basedocs.jl +++ b/base/docs/basedocs.jl @@ -281,6 +281,16 @@ kw"'" 1+1im 2+1im 3+1im 4+1im + julia> v = [1,2,3] + 3-element Array{Int64,1}: + 1 + 2 + 3 + + julia> v.' + 1×3 RowVector{Int64,Array{Int64,1}}: + 1 2 3 + """ kw".'" diff --git a/base/exports.jl b/base/exports.jl index 81caba8fb65c7..f14b89fbb6746 100644 --- a/base/exports.jl +++ b/base/exports.jl @@ -98,6 +98,7 @@ export RoundNearestTiesUp, RoundToZero, RoundUp, + RowVector, AbstractSerializer, SerializationState, Set, diff --git a/base/linalg/bidiag.jl b/base/linalg/bidiag.jl index 98db1a3fa26af..f040cd127c27d 100644 --- a/base/linalg/bidiag.jl +++ b/base/linalg/bidiag.jl @@ -354,6 +354,15 @@ A_mul_B!(C::AbstractVector, A::BiTri, B::AbstractVector) = A_mul_B_td!(C, A, B) A_mul_B!(C::AbstractMatrix, A::BiTri, B::AbstractVecOrMat) = A_mul_B_td!(C, A, B) A_mul_B!(C::AbstractVecOrMat, A::BiTri, B::AbstractVecOrMat) = A_mul_B_td!(C, A, B) +\(::Diagonal, ::RowVector) = throw(DimensionMismatch("Cannot left-divide matrix by transposed vector")) +\(::Bidiagonal, ::RowVector) = throw(DimensionMismatch("Cannot left-divide matrix by transposed vector")) +\{TA<:Number,TB<:Number}(::Bidiagonal{TA}, ::RowVector{TB}) = throw(DimensionMismatch("Cannot left-divide matrix by transposed vector")) + +At_ldiv_B(::Bidiagonal, ::RowVector) = throw(DimensionMismatch("Cannot left-divide matrix by transposed vector")) +At_ldiv_B{TA<:Number,TB<:Number}(::Bidiagonal{TA}, ::RowVector{TB}) = throw(DimensionMismatch("Cannot left-divide matrix by transposed vector")) + +Ac_ldiv_B(::Bidiagonal, ::RowVector) = throw(DimensionMismatch("Cannot left-divide matrix by transposed vector")) +Ac_ldiv_B{TA<:Number,TB<:Number}(::Bidiagonal{TA}, ::RowVector{TB}) = throw(DimensionMismatch("Cannot left-divide matrix by transposed vector")) function check_A_mul_B!_sizes(C, A, B) nA, mA = size(A) diff --git a/base/linalg/bitarray.jl b/base/linalg/bitarray.jl index 72124add2aac0..2b053e52745de 100644 --- a/base/linalg/bitarray.jl +++ b/base/linalg/bitarray.jl @@ -212,10 +212,6 @@ function findmin(a::BitArray) return m, mi end -## transpose ## - -transpose(B::BitVector) = reshape(copy(B), 1, length(B)) - # fast 8x8 bit transpose from Henry S. Warrens's "Hacker's Delight" # http://www.hackersdelight.org/hdcodetxt/transpose8.c.txt function transpose8x8(x::UInt64) diff --git a/base/linalg/dense.jl b/base/linalg/dense.jl index b4098fc2ddedf..e0e46de20c452 100644 --- a/base/linalg/dense.jl +++ b/base/linalg/dense.jl @@ -291,7 +291,7 @@ julia> kron(A, B) ``` """ function kron{T,S}(a::AbstractMatrix{T}, b::AbstractMatrix{S}) - R = Array{promote_type(T,S)}(size(a,1)*size(b,1), size(a,2)*size(b,2)) + R = Array{promote_op(*,T,S)}(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) aij = a[i,j] diff --git a/base/linalg/diagonal.jl b/base/linalg/diagonal.jl index cf81a2206a85b..cf56cdac76eee 100644 --- a/base/linalg/diagonal.jl +++ b/base/linalg/diagonal.jl @@ -254,6 +254,18 @@ function A_ldiv_B!{T}(D::Diagonal{T}, V::AbstractMatrix{T}) V end +# Methods to resolve ambiguities with `Diagonal` +@inline *(rowvec::RowVector, D::Diagonal) = transpose(D * transpose(rowvec)) +*(::Diagonal, ::RowVector) = throw(DimensionMismatch("Cannot right-multiply matrix by transposed vector")) + +@inline A_mul_Bt(D::Diagonal, rowvec::RowVector) = D*transpose(rowvec) + +At_mul_B(rowvec::RowVector, ::Diagonal) = throw(DimensionMismatch("Cannot left-multiply matrix by vector")) + +@inline A_mul_Bc(D::Diagonal, rowvec::RowVector) = D*ctranspose(rowvec) + +Ac_mul_B(rowvec::RowVector, ::Diagonal) = throw(DimensionMismatch("Cannot left-multiply matrix by vector")) + conj(D::Diagonal) = Diagonal(conj(D.diag)) transpose(D::Diagonal) = D ctranspose(D::Diagonal) = conj(D) diff --git a/base/linalg/generic.jl b/base/linalg/generic.jl index dff7113069929..9772fb7a0aa80 100644 --- a/base/linalg/generic.jl +++ b/base/linalg/generic.jl @@ -545,6 +545,20 @@ end @inline norm(x::Number, p::Real=2) = vecnorm(x, p) +@inline norm{T}(tv::RowVector{T}) = norm(transpose(tv)) + +""" + norm(rowvector, [q = 2]) + +Takes the q-norm of a `RowVector`, which is equivalent to the p-norm with +value `p = q/(q-1)`. They coincide at `p = q = 2`. + +The difference in norm between a vector space and its dual arises to preserve +the relationship between duality and the inner product, and the result is +consistent with the p-norm of `1 × n` matrix. +""" +@inline norm{T}(tv::RowVector{T}, q::Real) = q == Inf ? norm(transpose(tv), 1) : norm(transpose(tv), q/(q-1)) + function vecdot(x::AbstractArray, y::AbstractArray) lx = _length(x) if lx != _length(y) diff --git a/base/linalg/linalg.jl b/base/linalg/linalg.jl index e817b3a83d118..2c0d94b7c19e0 100644 --- a/base/linalg/linalg.jl +++ b/base/linalg/linalg.jl @@ -5,12 +5,13 @@ module LinAlg import Base: \, /, *, ^, +, -, == import Base: A_mul_Bt, At_ldiv_Bt, A_rdiv_Bc, At_ldiv_B, Ac_mul_Bc, A_mul_Bc, Ac_mul_B, Ac_ldiv_B, Ac_ldiv_Bc, At_mul_Bt, A_rdiv_Bt, At_mul_B -import Base: USE_BLAS64, abs, big, ceil, conj, convert, copy, copy!, ctranspose, - eltype, eye, findmax, findmin, fill!, floor, full, getindex, imag, inv, - isapprox, kron, ndims, parent, power_by_squaring, print_matrix, - promote_rule, real, round, setindex!, show, similar, size, transpose, trunc, - broadcast -using Base: promote_op, _length, iszero +import Base: USE_BLAS64, abs, big, broadcast, ceil, conj, convert, copy, copy!, + ctranspose, eltype, eye, findmax, findmin, fill!, floor, full, getindex, + hcat, imag, indices, inv, isapprox, kron, length, linearindexing, map, + ndims, parent, power_by_squaring, print_matrix, promote_rule, real, round, + setindex!, show, similar, size, transpose, trunc, typed_hcat +using Base: promote_op, _length, iszero, @pure, @propagate_inbounds, LinearFast, + reduce, hvcat_fill, typed_vcat, promote_typeof # 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))`. @@ -20,6 +21,7 @@ export BLAS, # Types + RowVector, SymTridiagonal, Tridiagonal, Bidiagonal, @@ -130,6 +132,7 @@ export trace, transpose, transpose!, + transpose_type, tril, triu, tril!, @@ -235,6 +238,7 @@ copy_oftype{T,N}(A::AbstractArray{T,N}, ::Type{T}) = copy(A) copy_oftype{T,N,S}(A::AbstractArray{T,N}, ::Type{S}) = convert(AbstractArray{S,N}, A) include("transpose.jl") +include("rowvector.jl") include("exceptions.jl") include("generic.jl") diff --git a/base/linalg/matmul.jl b/base/linalg/matmul.jl index 43c49d9c1d250..4dbd4c8dae817 100644 --- a/base/linalg/matmul.jl +++ b/base/linalg/matmul.jl @@ -69,9 +69,7 @@ function dot{T<:BlasComplex, TI<:Integer}(x::Vector{T}, rx::Union{UnitRange{TI}, BLAS.dotc(length(rx), pointer(x)+(first(rx)-1)*sizeof(T), step(rx), pointer(y)+(first(ry)-1)*sizeof(T), step(ry)) end -Ac_mul_B(x::AbstractVector, y::AbstractVector) = [dot(x, y)] -At_mul_B{T<:Real}(x::AbstractVector{T}, y::AbstractVector{T}) = [dot(x, y)] -At_mul_B{T<:BlasComplex}(x::StridedVector{T}, y::StridedVector{T}) = [BLAS.dotu(x, y)] +At_mul_B{T<:BlasComplex}(x::StridedVector{T}, y::StridedVector{T}) = BLAS.dotu(x, y) # Matrix-vector multiplication function (*){T<:BlasFloat,S}(A::StridedMatrix{T}, x::StridedVector{S}) @@ -82,7 +80,6 @@ function (*){T,S}(A::AbstractMatrix{T}, x::AbstractVector{S}) TS = promote_op(matprod, T, S) A_mul_B!(similar(x,TS,size(A,1)),A,x) end -(*)(A::AbstractVector, B::AbstractMatrix) = reshape(A,length(A),1)*B A_mul_B!{T<:BlasFloat}(y::StridedVector{T}, A::StridedVecOrMat{T}, x::StridedVector{T}) = gemv!(y, 'N', A, x) for elty in (Float32,Float64) diff --git a/base/linalg/qr.jl b/base/linalg/qr.jl index d3acc6a562b1f..0e53bf6572219 100644 --- a/base/linalg/qr.jl +++ b/base/linalg/qr.jl @@ -543,6 +543,8 @@ function A_mul_Bc(A::AbstractMatrix, B::Union{QRCompactWYQ,QRPackedQ}) throw(DimensionMismatch("matrix A has dimensions $(size(A)) but matrix B has dimensions $(size(B))")) end end +@inline A_mul_Bc(rowvec::RowVector, B::Union{LinAlg.QRCompactWYQ,LinAlg.QRPackedQ}) = ctranspose(B*ctranspose(rowvec)) + ### AcQ/AcQc for (f1, f2) in ((:Ac_mul_B, :A_mul_B!), diff --git a/base/linalg/rowvector.jl b/base/linalg/rowvector.jl new file mode 100644 index 0000000000000..a87adf0be459a --- /dev/null +++ b/base/linalg/rowvector.jl @@ -0,0 +1,214 @@ +""" + RowVector(vector) + +A lazy-view wrapper of an `AbstractVector`, which turns a length-`n` vector into +a `1×n` shaped row vector and represents the transpose of a vector (the elements +are also transposed recursively). This type is usually constructed (and +unwrapped) via the `transpose()` function or `.'` operator (or related +`ctranspose()` or `'` operator). + +By convention, a vector can be multiplied by a matrix on its left (`A * v`) +whereas a row vector can be multiplied by a matrix on its right (such that +`v.' * A = (A.' * v).'`). It differs from a `1×n`-sized matrix by the facts that +its transpose returns a vector and the inner product `v1.' * v2` returns a +scalar, but will otherwise behave similarly. +""" +immutable RowVector{T,V<:AbstractVector} <: AbstractMatrix{T} + vec::V + function RowVector(v::V) + check_types(T,v) + new(v) + end +end + + +@inline check_types{T1,T2}(::Type{T1},::AbstractVector{T2}) = check_types(T1, T2) +@pure check_types{T1,T2}(::Type{T1},::Type{T2}) = T1 === transpose_type(T2) ? nothing : error("Element type mismatch. Tried to create a `RowVector{$T1}` from an `AbstractVector{$T2}`") + +# The element type may be transformed as transpose is recursive +@inline transpose_type{T}(::Type{T}) = promote_op(transpose, T) + +# Constructors that take a vector +@inline RowVector{T}(vec::AbstractVector{T}) = RowVector{transpose_type(T),typeof(vec)}(vec) +@inline (::Type{RowVector{T}}){T}(vec::AbstractVector{T}) = RowVector{T,typeof(vec)}(vec) + +# Constructors that take a size and default to Array +@inline (::Type{RowVector{T}}){T}(n::Int) = RowVector{T}(Vector{transpose_type(T)}(n)) +@inline (::Type{RowVector{T}}){T}(n1::Int, n2::Int) = n1 == 1 ? RowVector{T}(Vector{transpose_type(T)}(n2)) : error("RowVector expects 1×N size, got ($n1,$n2)") +@inline (::Type{RowVector{T}}){T}(n::Tuple{Int}) = RowVector{T}(Vector{transpose_type(T)}(n[1])) +@inline (::Type{RowVector{T}}){T}(n::Tuple{Int,Int}) = n[1] == 1 ? RowVector{T}(Vector{transpose_type(T)}(n[2])) : error("RowVector expects 1×N size, got $n") + +# Conversion of underlying storage +convert{T,V<:AbstractVector}(::Type{RowVector{T,V}}, rowvec::RowVector) = RowVector{T,V}(convert(V,rowvec.vec)) + +# similar() +@inline similar(rowvec::RowVector) = RowVector(similar(rowvec.vec)) +@inline similar{T}(rowvec::RowVector, ::Type{T}) = RowVector(similar(rowvec.vec, transpose_type(T))) +# There is no resizing similar() because it would be ambiguous if the result were a Matrix or a RowVector + +# Basic methods +""" + transpose(v::AbstractVector) + +The transposition operator (`.'`). + +# Example + +```jldoctest +julia> v = [1,2,3] +3-element Array{Int64,1}: + 1 + 2 + 3 + +julia> transpose(v) +1×3 RowVector{Int64,Array{Int64,1}}: + 1 2 3 +``` +""" +@inline transpose(vec::AbstractVector) = RowVector(vec) +@inline ctranspose{T}(vec::AbstractVector{T}) = RowVector(conj(vec)) +@inline ctranspose{T<:Real}(vec::AbstractVector{T}) = RowVector(vec) + +@inline transpose(rowvec::RowVector) = rowvec.vec +@inline ctranspose{T}(rowvec::RowVector{T}) = conj(rowvec.vec) +@inline ctranspose{T<:Real}(rowvec::RowVector{T}) = rowvec.vec + +parent(rowvec::RowVector) = rowvec.vec + +# Strictly, these are unnecessary but will make things stabler if we introduce +# a "view" for conj(::AbstractArray) +@inline conj(rowvec::RowVector) = RowVector(conj(rowvec.vec)) +@inline conj{T<:Real}(rowvec::RowVector{T}) = rowvec + +# AbstractArray interface +@inline length(rowvec::RowVector) = length(rowvec.vec) +@inline size(rowvec::RowVector) = (1, length(rowvec.vec)) +@inline size(rowvec::RowVector, d) = ifelse(d==2, length(rowvec.vec), 1) +@inline indices(rowvec::RowVector) = (Base.OneTo(1), indices(rowvec.vec)[1]) +@inline indices(rowvec::RowVector, d) = ifelse(d == 2, indices(rowvec.vec)[1], Base.OneTo(1)) +linearindexing{V<:RowVector}(::Union{V,Type{V}}) = LinearFast() + +@propagate_inbounds getindex(rowvec::RowVector, i) = transpose(rowvec.vec[i]) +@propagate_inbounds setindex!(rowvec::RowVector, v, i) = setindex!(rowvec.vec, transpose(v), i) + +# Cartesian indexing is distorted by getindex +# Furthermore, Cartesian indexes don't have to match shape, apparently! +@inline function getindex(rowvec::RowVector, i::CartesianIndex) + @boundscheck if !(i.I[1] == 1 && i.I[2] ∈ indices(rowvec.vec)[1] && check_tail_indices(i.I...)) + throw(BoundsError(rowvec, i.I)) + end + @inbounds return transpose(rowvec.vec[i.I[2]]) +end +@inline function setindex!(rowvec::RowVector, v, i::CartesianIndex) + @boundscheck if !(i.I[1] == 1 && i.I[2] ∈ indices(rowvec.vec)[1] && check_tail_indices(i.I...)) + throw(BoundsError(rowvec, i.I)) + end + @inbounds rowvec.vec[i.I[2]] = transpose(v) +end + +@propagate_inbounds getindex(rowvec::RowVector, ::CartesianIndex{0}) = getindex(rowvec) +@propagate_inbounds getindex(rowvec::RowVector, i::CartesianIndex{1}) = getindex(rowvec, i.I[1]) + +@propagate_inbounds setindex!(rowvec::RowVector, v, ::CartesianIndex{0}) = setindex!(rowvec, v) +@propagate_inbounds setindex!(rowvec::RowVector, v, i::CartesianIndex{1}) = setindex!(rowvec, v, i.I[1]) + +@inline check_tail_indices(i1, i2) = true +@inline check_tail_indices(i1, i2, i3, is...) = i3 == 1 ? check_tail_indices(i1, i2, is...) : false + +# helper function for below +@inline to_vec(rowvec::RowVector) = transpose(rowvec) +@inline to_vec(x::Number) = x +@inline to_vecs(rowvecs...) = (map(to_vec, rowvecs)...) + +# map +@inline map(f, rowvecs::RowVector...) = RowVector(map(f, to_vecs(rowvecs...)...)) + +# broacast (other combinations default to higher-dimensional array) +@inline broadcast(f, rowvecs::Union{Number,RowVector}...) = RowVector(broadcast(f, to_vecs(rowvecs...)...)) + +# Horizontal concatenation # + +@inline hcat(X::RowVector...) = transpose(vcat(map(transpose, X)...)) +@inline hcat(X::Union{RowVector,Number}...) = transpose(vcat(map(transpose, X)...)) + +@inline typed_hcat{T}(::Type{T}, X::RowVector...) = transpose(typed_vcat(T, map(transpose, X)...)) +@inline typed_hcat{T}(::Type{T}, X::Union{RowVector,Number}...) = transpose(typed_vcat(T, map(transpose, X)...)) + +# Multiplication # + +@inline function *(rowvec::RowVector, vec::AbstractVector) + if length(rowvec) != length(vec) + throw(DimensionMismatch("A has dimensions $(size(rowvec)) but B has dimensions $(size(vec))")) + end + sum(@inbounds(return rowvec[i]*vec[i]) for i = 1:length(vec)) +end +@inline *(rowvec::RowVector, mat::AbstractMatrix) = transpose(mat.' * transpose(rowvec)) +*(vec::AbstractVector, mat::AbstractMatrix) = throw(DimensionMismatch("Cannot left-multiply a matrix by a vector")) # Should become a deprecation +*(::RowVector, ::RowVector) = throw(DimensionMismatch("Cannot multiply two transposed vectors")) +@inline *(vec::AbstractVector, rowvec::RowVector) = vec .* rowvec +*(vec::AbstractVector, rowvec::AbstractVector) = throw(DimensionMismatch("Cannot multiply two vectors")) +*(mat::AbstractMatrix, rowvec::RowVector) = throw(DimensionMismatch("Cannot right-multiply matrix by transposed vector")) + +# Transposed forms +A_mul_Bt(::RowVector, ::AbstractVector) = throw(DimensionMismatch("Cannot multiply two transposed vectors")) +@inline A_mul_Bt(rowvec::RowVector, mat::AbstractMatrix) = transpose(mat * transpose(rowvec)) +A_mul_Bt(vec::AbstractVector, mat::AbstractMatrix) = throw(DimensionMismatch("Cannot left-multiply a matrix by a vector")) +@inline A_mul_Bt(rowvec1::RowVector, rowvec2::RowVector) = rowvec1*transpose(rowvec2) +A_mul_Bt(vec::AbstractVector, rowvec::RowVector) = throw(DimensionMismatch("Cannot multiply two vectors")) +@inline A_mul_Bt(vec1::AbstractVector, vec2::AbstractVector) = vec1 * transpose(vec2) +@inline A_mul_Bt(mat::AbstractMatrix, rowvec::RowVector) = mat * transpose(rowvec) + +@inline At_mul_Bt(rowvec::RowVector, vec::AbstractVector) = transpose(rowvec) * transpose(vec) +At_mul_Bt(rowvec::RowVector, mat::AbstractMatrix) = throw(DimensionMismatch("Cannot left-multiply matrix by vector")) +@inline At_mul_Bt(vec::AbstractVector, mat::AbstractMatrix) = transpose(mat * vec) +At_mul_Bt(rowvec1::RowVector, rowvec2::RowVector) = throw(DimensionMismatch("Cannot multiply two vectors")) +@inline At_mul_Bt(vec::AbstractVector, rowvec::RowVector) = transpose(vec)*transpose(rowvec) +At_mul_Bt(vec::AbstractVector, rowvec::AbstractVector) = throw(DimensionMismatch("Cannot multiply two transposed vectors")) +@inline At_mul_Bt(mat::AbstractMatrix, rowvec::RowVector) = mat.' * transpose(rowvec) + +At_mul_B(::RowVector, ::AbstractVector) = throw(DimensionMismatch("Cannot multiply two vectors")) +At_mul_B(rowvec::RowVector, mat::AbstractMatrix) = throw(DimensionMismatch("Cannot left-multiply matrix by vector")) +@inline At_mul_B(vec::AbstractVector, mat::AbstractMatrix) = transpose(At_mul_B(mat,vec)) +@inline At_mul_B(rowvec1::RowVector, rowvec2::RowVector) = transpose(rowvec1) * rowvec2 +At_mul_B(vec::AbstractVector, rowvec::RowVector) = throw(DimensionMismatch("Cannot multiply two transposed vectors")) +@inline At_mul_B{T<:Real}(vec1::AbstractVector{T}, vec2::AbstractVector{T}) = reduce(+, map(At_mul_B, vec1, vec2)) # Seems to be overloaded... +@inline At_mul_B(vec1::AbstractVector, vec2::AbstractVector) = transpose(vec1) * vec2 +At_mul_B(mat::AbstractMatrix, rowvec::RowVector) = throw(DimensionMismatch("Cannot right-multiply matrix by transposed vector")) + +# Conjugated forms +A_mul_Bc(::RowVector, ::AbstractVector) = throw(DimensionMismatch("Cannot multiply two transposed vectors")) +@inline A_mul_Bc(rowvec::RowVector, mat::AbstractMatrix) = ctranspose(mat * ctranspose(rowvec)) +A_mul_Bc(vec::AbstractVector, mat::AbstractMatrix) = throw(DimensionMismatch("Cannot left-multiply a matrix by a vector")) +@inline A_mul_Bc(rowvec1::RowVector, rowvec2::RowVector) = rowvec1 * ctranspose(rowvec2) +A_mul_Bc(vec::AbstractVector, rowvec::RowVector) = throw(DimensionMismatch("Cannot multiply two vectors")) +@inline A_mul_Bc(vec1::AbstractVector, vec2::AbstractVector) = vec1 * ctranspose(vec2) +@inline A_mul_Bc(mat::AbstractMatrix, rowvec::RowVector) = mat * ctranspose(rowvec) + +@inline Ac_mul_Bc(rowvec::RowVector, vec::AbstractVector) = ctranspose(rowvec) * ctranspose(vec) +Ac_mul_Bc(rowvec::RowVector, mat::AbstractMatrix) = throw(DimensionMismatch("Cannot left-multiply matrix by vector")) +@inline Ac_mul_Bc(vec::AbstractVector, mat::AbstractMatrix) = ctranspose(mat * vec) +Ac_mul_Bc(rowvec1::RowVector, rowvec2::RowVector) = throw(DimensionMismatch("Cannot multiply two vectors")) +@inline Ac_mul_Bc(vec::AbstractVector, rowvec::RowVector) = ctranspose(vec)*ctranspose(rowvec) +Ac_mul_Bc(vec::AbstractVector, rowvec::AbstractVector) = throw(DimensionMismatch("Cannot multiply two transposed vectors")) +@inline Ac_mul_Bc(mat::AbstractMatrix, rowvec::RowVector) = mat' * ctranspose(rowvec) + +Ac_mul_B(::RowVector, ::AbstractVector) = throw(DimensionMismatch("Cannot multiply two vectors")) +Ac_mul_B(rowvec::RowVector, mat::AbstractMatrix) = throw(DimensionMismatch("Cannot left-multiply matrix by vector")) +@inline Ac_mul_B(vec::AbstractVector, mat::AbstractMatrix) = ctranspose(Ac_mul_B(mat,vec)) +@inline Ac_mul_B(rowvec1::RowVector, rowvec2::RowVector) = ctranspose(rowvec1) * rowvec2 +Ac_mul_B(vec::AbstractVector, rowvec::RowVector) = throw(DimensionMismatch("Cannot multiply two transposed vectors")) +@inline Ac_mul_B(vec1::AbstractVector, vec2::AbstractVector) = ctranspose(vec1)*vec2 +Ac_mul_B(mat::AbstractMatrix, rowvec::RowVector) = throw(DimensionMismatch("Cannot right-multiply matrix by transposed vector")) + +# Left Division # + +\(mat::AbstractMatrix, rowvec::RowVector) = throw(DimensionMismatch("Cannot left-divide transposed vector by matrix")) +At_ldiv_B(mat::AbstractMatrix, rowvec::RowVector) = throw(DimensionMismatch("Cannot left-divide transposed vector by matrix")) +Ac_ldiv_B(mat::AbstractMatrix, rowvec::RowVector) = throw(DimensionMismatch("Cannot left-divide transposed vector by matrix")) + +# Right Division # + +@inline /(rowvec::RowVector, mat::AbstractMatrix) = transpose(transpose(mat) \ transpose(rowvec)) +@inline A_rdiv_Bt(rowvec::RowVector, mat::AbstractMatrix) = transpose(mat \ transpose(rowvec)) +@inline A_rdiv_Bc(rowvec::RowVector, mat::AbstractMatrix) = ctranspose(mat \ ctranspose(rowvec)) diff --git a/base/linalg/transpose.jl b/base/linalg/transpose.jl index ab4007486509b..7210c3198dc45 100644 --- a/base/linalg/transpose.jl +++ b/base/linalg/transpose.jl @@ -128,9 +128,6 @@ end @inline ctranspose{T<:Real}(A::AbstractVecOrMat{T}) = transpose(A) -transpose(x::AbstractVector) = [ transpose(v) for i=Base.of_indices(x, Base.OneTo(1)), v in x ] -ctranspose{T}(x::AbstractVector{T}) = T[ ctranspose(v) for i=Base.of_indices(x, Base.OneTo(1)), v in x ] - function copy_transpose!{R,S}(B::AbstractVecOrMat{R}, ir_dest::Range{Int}, jr_dest::Range{Int}, A::AbstractVecOrMat{S}, ir_src::Range{Int}, jr_src::Range{Int}) if length(ir_dest) != length(jr_src) diff --git a/base/linalg/triangular.jl b/base/linalg/triangular.jl index a0eda4499fdf6..3af8ab2dc2b85 100644 --- a/base/linalg/triangular.jl +++ b/base/linalg/triangular.jl @@ -1634,6 +1634,47 @@ At_mul_Bt(A::AbstractTriangular, B::AbstractTriangular) = At_mul_B(A, B.') At_mul_Bt(A::AbstractTriangular, B::AbstractMatrix) = At_mul_B(A, B.') At_mul_Bt(A::AbstractMatrix, B::AbstractTriangular) = A_mul_Bt(A.', B) +# Specializations for RowVector +@inline *(rowvec::RowVector, A::AbstractTriangular) = transpose(A * transpose(rowvec)) +*(::AbstractTriangular, ::RowVector) = throw(DimensionMismatch("Cannot right-multiply matrix by transposed vector")) + +@inline A_mul_Bt(rowvec::RowVector, A::AbstractTriangular) = transpose(A * transpose(rowvec)) +@inline A_mul_Bt(A::AbstractTriangular, rowvec::RowVector) = A * transpose(rowvec) + +@inline At_mul_Bt(A::AbstractTriangular, rowvec::RowVector) = A.' * transpose(rowvec) +At_mul_Bt(::RowVector, ::AbstractTriangular) = throw(DimensionMismatch("Cannot left-multiply matrix by vector")) + +At_mul_B(::AbstractTriangular, ::RowVector) = throw(DimensionMismatch("Cannot right-multiply matrix by transposed vector")) +At_mul_B(::RowVector, ::AbstractTriangular) = throw(DimensionMismatch("Cannot left-multiply matrix by vector")) + +@inline A_mul_Bc(rowvec::RowVector, A::AbstractTriangular) = ctranspose(A * ctranspose(rowvec)) +@inline A_mul_Bc(A::AbstractTriangular, rowvec::RowVector) = A * ctranspose(rowvec) + +@inline Ac_mul_Bc(A::AbstractTriangular, rowvec::RowVector) = A' * ctranspose(rowvec) +Ac_mul_Bc(::RowVector, ::AbstractTriangular) = throw(DimensionMismatch("Cannot left-multiply matrix by vector")) + +Ac_mul_B(::AbstractTriangular, ::RowVector) = throw(DimensionMismatch("Cannot right-multiply matrix by transposed vector")) +Ac_mul_B(::RowVector, ::AbstractTriangular) = throw(DimensionMismatch("Cannot left-multiply matrix by vector")) + +@inline /(rowvec::RowVector, A::Union{UpperTriangular,LowerTriangular}) = transpose(transpose(A) \ transpose(rowvec)) +@inline /(rowvec::RowVector, A::Union{UnitUpperTriangular,UnitLowerTriangular}) = transpose(transpose(A) \ transpose(rowvec)) + +@inline A_rdiv_Bt(rowvec::RowVector, A::Union{UpperTriangular,LowerTriangular}) = transpose(A \ transpose(rowvec)) +@inline A_rdiv_Bt(rowvec::RowVector, A::Union{UnitUpperTriangular,UnitLowerTriangular}) = transpose(A \ transpose(rowvec)) + +@inline A_rdiv_Bc(rowvec::RowVector, A::Union{UpperTriangular,LowerTriangular}) = ctranspose(A \ ctranspose(rowvec)) +@inline A_rdiv_Bc(rowvec::RowVector, A::Union{UnitUpperTriangular,UnitLowerTriangular}) = ctranspose(A \ ctranspose(rowvec)) + +\(::Union{UpperTriangular,LowerTriangular}, ::RowVector) = throw(DimensionMismatch("Cannot left-divide matrix by transposed vector")) +\(::Union{UnitUpperTriangular,UnitLowerTriangular}, ::RowVector) = throw(DimensionMismatch("Cannot left-divide matrix by transposed vector")) + +At_ldiv_B(::Union{UpperTriangular,LowerTriangular}, ::RowVector) = throw(DimensionMismatch("Cannot left-divide matrix by transposed vector")) +At_ldiv_B(::Union{UnitUpperTriangular,UnitLowerTriangular}, ::RowVector) = throw(DimensionMismatch("Cannot left-divide matrix by transposed vector")) + +Ac_ldiv_B(::Union{UpperTriangular,LowerTriangular}, ::RowVector) = throw(DimensionMismatch("Cannot left-divide matrix by transposed vector")) +Ac_ldiv_B(::Union{UnitUpperTriangular,UnitLowerTriangular}, ::RowVector) = throw(DimensionMismatch("Cannot left-divide matrix by transposed vector")) + + # Complex matrix logarithm for the upper triangular factor, see: # Al-Mohy and Higham, "Improved inverse scaling and squaring algorithms for # the matrix logarithm", SIAM J. Sci. Comput., 34(4), (2012), pp. C153-C169. diff --git a/base/sparse/cholmod.jl b/base/sparse/cholmod.jl index 5bccc8f270c26..3cf3c0fc4d331 100644 --- a/base/sparse/cholmod.jl +++ b/base/sparse/cholmod.jl @@ -1570,6 +1570,7 @@ function (\)(L::FactorComponent, B::SparseVecOrMat) end Ac_ldiv_B(L::FactorComponent, B) = ctranspose(L)\B +Ac_ldiv_B(L::FactorComponent, B::RowVector) = ctranspose(L)\B # ambiguity (\){T<:VTypes}(L::Factor{T}, B::Dense{T}) = solve(CHOLMOD_A, L, B) (\)(L::Factor{Float64}, B::VecOrMat{Complex{Float64}}) = complex.(L\real(B), L\imag(B)) diff --git a/base/sparse/linalg.jl b/base/sparse/linalg.jl index 18570640c655d..1345c9fa848e9 100644 --- a/base/sparse/linalg.jl +++ b/base/sparse/linalg.jl @@ -872,6 +872,8 @@ function (\)(A::SparseMatrixCSC, B::AbstractVecOrMat) end end +(\)(::SparseMatrixCSC, ::RowVector) = throw(DimensionMismatch("Cannot left-divide matrix by transposed vector")) + function factorize(A::SparseMatrixCSC) m, n = size(A) if m == n diff --git a/base/sparse/sparsevector.jl b/base/sparse/sparsevector.jl index 2fda783c40b9c..63c40f7d79a7e 100644 --- a/base/sparse/sparsevector.jl +++ b/base/sparse/sparsevector.jl @@ -1316,18 +1316,8 @@ vecnorm(x::AbstractSparseVector, p::Real=2) = vecnorm(nonzeros(x), p) # Transpose # (The only sparse matrix structure in base is CSC, so a one-row sparse matrix is worse than dense) -transpose(x::SparseVector) = _ct(identity, x) -ctranspose(x::SparseVector) = _ct(conj, x) -function _ct{T}(f, x::SparseVector{T}) - isempty(x) && return Array{T}(1, 0) - A = zeros(T, 1, length(x)) - xnzind = nonzeroinds(x) - xnzval = nonzeros(x) - for (i,v) in zip(xnzind, xnzval) - @inbounds A[i] = f(v) - end - A -end +@inline transpose(sv::SparseVector) = RowVector(sv) +@inline ctranspose(sv::SparseVector) = RowVector(conj(sv)) ### BLAS Level-1 diff --git a/doc/src/stdlib/linalg.md b/doc/src/stdlib/linalg.md index 06808492c42bf..a3a3a74cb8bd9 100644 --- a/doc/src/stdlib/linalg.md +++ b/doc/src/stdlib/linalg.md @@ -99,6 +99,7 @@ Base.LinAlg.istril Base.LinAlg.istriu Base.LinAlg.isdiag Base.LinAlg.ishermitian +Base.LinAlg.RowVector Base.transpose Base.transpose! Base.ctranspose diff --git a/test/arrayops.jl b/test/arrayops.jl index b2cd9fb77b861..a46d640b26a5a 100644 --- a/test/arrayops.jl +++ b/test/arrayops.jl @@ -591,13 +591,13 @@ end @test R == [1, 1, 2, 2, 1, 1, 2, 2] R = repeat([1, 2], inner = (1, 1), outer = (1, 1)) - @test R == [1, 2]'' + @test R == reshape([1, 2], (2,1)) R = repeat([1, 2], inner = (2, 1), outer = (1, 1)) - @test R == [1, 1, 2, 2]'' + @test R == reshape([1, 1, 2, 2], (4,1)) R = repeat([1, 2], inner = (1, 2), outer = (1, 1)) @test R == [1 1; 2 2] R = repeat([1, 2], inner = (1, 1), outer = (2, 1)) - @test R == [1, 2, 1, 2]'' + @test R == reshape([1, 2, 1, 2], (4,1)) R = repeat([1, 2], inner = (1, 1), outer = (1, 2)) @test R == [1 1; 2 2] diff --git a/test/bitarray.jl b/test/bitarray.jl index 8413cd5bfc39e..8e88c2dd8aef5 100644 --- a/test/bitarray.jl +++ b/test/bitarray.jl @@ -7,6 +7,7 @@ using Base: findprevnot, findnextnot tc{N}(r1::NTuple{N}, r2::NTuple{N}) = all(x->tc(x...), [zip(r1,r2)...]) tc{N}(r1::BitArray{N}, r2::Union{BitArray{N},Array{Bool,N}}) = true +tc(r1::RowVector{Bool,BitVector}, r2::Union{RowVector{Bool,BitVector},RowVector{Bool,Vector{Bool}}}) = true tc{T}(r1::T, r2::T) = true tc(r1,r2) = false @@ -781,7 +782,7 @@ timesofar("unary arithmetic") @check_bit_operation broadcast(|, b0, b0) BitVector @check_bit_operation broadcast(xor, b0, b0) BitVector @check_bit_operation broadcast(*, b0, b0) BitVector - @check_bit_operation (*)(b0, b0') Matrix{Int} + @check_bit_operation (*)(b0, b0') BitMatrix end @testset "Matrix{Bool}/Matrix{Int}" begin @@ -1284,7 +1285,7 @@ end @testset "Transpose" begin b1 = bitrand(v1) - @check_bit_operation transpose(b1) BitMatrix + @check_bit_operation transpose(b1) RowVector{Bool,BitVector} for m1 = 0:n1, m2 = 0:n2 b1 = bitrand(m1, m2) diff --git a/test/choosetests.jl b/test/choosetests.jl index cb1f045e72aa1..55689a9121a54 100644 --- a/test/choosetests.jl +++ b/test/choosetests.jl @@ -130,7 +130,7 @@ function choosetests(choices = []) "linalg/diagonal", "linalg/pinv", "linalg/givens", "linalg/cholesky", "linalg/lu", "linalg/symmetric", "linalg/generic", "linalg/uniformscaling", "linalg/lq", - "linalg/hessenberg"] + "linalg/hessenberg", "linalg/rowvector"] if Base.USE_GPL_LIBS push!(linalgtests, "linalg/arnoldi") end diff --git a/test/linalg/dense.jl b/test/linalg/dense.jl index f3c14beb9b779..28661f7200ee4 100644 --- a/test/linalg/dense.jl +++ b/test/linalg/dense.jl @@ -392,7 +392,7 @@ let end # similar issue for Array{Real} -@test Real[1 2] * Real[1.5; 2.0] == [5.5] +@test Real[1 2] * Real[1.5; 2.0] == Real[5.5] # Matrix exponential for elty in (Float32, Float64, Complex64, Complex128) diff --git a/test/linalg/matmul.jl b/test/linalg/matmul.jl index 79810e8de329b..44421f3190e92 100644 --- a/test/linalg/matmul.jl +++ b/test/linalg/matmul.jl @@ -256,7 +256,7 @@ for elty in (Float32,Float64,Complex64,Complex128) @test_throws BoundsError dot(x, 1:4, y, 1:4) @test_throws BoundsError dot(x, 1:3, y, 2:4) @test dot(x,1:2,y,1:2) == convert(elty,12.5) - @test x.'*y == convert(Vector{elty},[29.0]) + @test x.'*y == convert(elty,29.0) end vecdot_(x,y) = invoke(vecdot, Tuple{Any,Any}, x,y) # generic vecdot diff --git a/test/linalg/rowvector.jl b/test/linalg/rowvector.jl new file mode 100644 index 0000000000000..8b23a8ed7ebdc --- /dev/null +++ b/test/linalg/rowvector.jl @@ -0,0 +1,245 @@ +# This file is a part of Julia. License is MIT: http://julialang.org/license + +@testset "RowVector" begin + +@testset "Core" begin + v = [1,2,3] + z = [1+im,2,3] + + @test RowVector(v) == [1 2 3] + @test RowVector{Int}(v) == [1 2 3] + @test size(RowVector{Int}(3)) === (1,3) + @test size(RowVector{Int}(1,3)) === (1,3) + @test size(RowVector{Int}((3,))) === (1,3) + @test size(RowVector{Int}((1,3))) === (1,3) + @test_throws Exception RowVector{Float64, Vector{Int}}(v) + + @test (v.')::RowVector == [1 2 3] + @test (v')::RowVector == [1 2 3] + @test (z.')::RowVector == [1+im 2 3] + @test (z')::RowVector == [1-im 2 3] + + rv = v.' + tz = z.' + + @test (rv.')::Vector == [1, 2, 3] + @test (rv')::Vector == [1, 2, 3] + @test (tz.')::Vector == [1+im, 2, 3] + @test (tz')::Vector == [1-im, 2, 3] + + @test conj(rv) === rv + @test conj(tz) == [1-im 2 3] + + @test isa(similar(rv), RowVector) + @test isa(similar(rv, Float64), RowVector) + @test isa(copy(rv), RowVector) + + @test rv[2] === v[2] + @test rv[1,2] === v[2] + + @test (rv2 = copy(rv); rv2[2] = 6; rv2[2] === 6) + @test (rv2 = copy(rv); rv2[1,2] = 6; rv2[2] === 6) + + @test length(rv) === 3 + @test size(rv) === (1,3) + @test size(rv,1) === 1 + @test size(rv,2) === 3 + + @test map(-, rv)::RowVector == [-1 -2 -3] + @test (-).(rv)::RowVector == [-1 -2 -3] + @test (-).(rv,1)::RowVector == [0 1 2] + + y = rand(Complex{Float64},3) + @test sum(abs2, imag.(diag(y .+ y'))) < 1e-20 +end + +@testset "Diagonal ambiguity methods" begin + d = Diagonal([1,2,3]) + v = [2,3,4] + rv = v.' + + @test (rv*d)::RowVector == [2,6,12].' + @test_throws Exception d*rv + + @test (d*rv.')::Vector == [2,6,12] + + @test_throws Exception rv.'*d + + @test (d*rv')::Vector == [2,6,12] + + @test_throws Exception rv'*d + + @test (rv/d)::RowVector ≈ [2/1 3/2 4/3] + + @test_throws Exception d \ rv + +end + +@testset "Bidiagonal ambiguity methods" begin + bd = Bidiagonal([1,2,3], [0,0], true) + v = [2,3,4] + rv = v.' + + @test (rv/bd)::RowVector ≈ [2/1 3/2 4/3] + + @test_throws Exception bd \ rv +end + +@testset "hcat" begin + @test isa([([1, 2, 3].') 4], RowVector{Int}) + @test isa([([1, 2, 3].') ([4, 5].')], RowVector{Int}) +end + +@testset "Left Division" begin + mat = diagm([1,2,3]) + v = [2,3,4] + rv = v.' + + @test_throws Exception mat \ rv +end + +@testset "Multiplication" begin + v = [1,2,3] + rv = v.' + mat = diagm([1,2,3]) + + @test (rv*v) === 14 + @test (rv*mat)::RowVector == [1 4 9] + @test_throws Exception [1]*reshape([1],(1,1)) # no longer permitted + @test_throws Exception rv*rv + @test (v*rv)::Matrix == [1 2 3; 2 4 6; 3 6 9] + @test_throws Exception v*v # Was previously a missing method error, now an error message + @test_throws Exception mat*rv + + @test_throws Exception rv*v.' + @test (rv*mat.')::RowVector == [1 4 9] + @test_throws Exception [1]*reshape([1],(1,1)).' # no longer permitted + @test rv*rv.' === 14 + @test_throws Exception v*rv.' + @test (v*v.')::Matrix == [1 2 3; 2 4 6; 3 6 9] + @test (mat*rv.')::Vector == [1,4,9] + + @test (rv.'*v.')::Matrix == [1 2 3; 2 4 6; 3 6 9] + @test_throws Exception rv.'*mat.' + @test (v.'*mat.')::RowVector == [1 4 9] + @test_throws Exception rv.'*rv.' + @test v.'*rv.' === 14 + @test_throws Exception v.'*v.' + @test (mat.'*rv.')::Vector == [1,4,9] + + @test_throws Exception rv.'*v + @test_throws Exception rv.'*mat + @test (v.'*mat)::RowVector == [1 4 9] + @test (rv.'*rv)::Matrix == [1 2 3; 2 4 6; 3 6 9] + @test_throws Exception v.'*rv + @test v.'*v === 14 + @test_throws Exception mat.'*rv + + z = [1+im,2,3] + cz = z' + mat = diagm([1+im,2,3]) + + @test cz*z === 15 + 0im + + @test_throws Exception cz*z' + @test (cz*mat')::RowVector == [-2im 4 9] + @test_throws Exception [1]*reshape([1],(1,1))' # no longer permitted + @test cz*cz' === 15 + 0im + @test_throws Exception z*vz' + @test (z*z')::Matrix == [2 2+2im 3+3im; 2-2im 4 6; 3-3im 6 9] + @test (mat*cz')::Vector == [2im,4,9] + + @test (cz'*z')::Matrix == [2 2+2im 3+3im; 2-2im 4 6; 3-3im 6 9] + @test_throws Exception cz'*mat' + @test (z'*mat')::RowVector == [-2im 4 9] + @test_throws Exception cz'*cz' + @test z'*cz' === 15 + 0im + @test_throws Exception z'*z' + @test (mat'*cz')::Vector == [2,4,9] + + @test_throws Exception cz'*z + @test_throws Exception cz'*mat + @test (z'*mat)::RowVector == [2 4 9] + @test (cz'*cz)::Matrix == [2 2+2im 3+3im; 2-2im 4 6; 3-3im 6 9] + @test_throws Exception z'*cz + @test z'*z === 15 + 0im + @test_throws Exception mat'*cz + +end + +@testset "norm" begin + @test norm([3.0,4.0].') ≈ 5.0 + @test norm([3.0,4.0].', 1) ≈ 4.0 + @test norm([3.0,4.0].', Inf) ≈ 7.0 +end + +@testset "QR ambiguity methods" begin + qrmat = Base.LinAlg.getq(qrfact(eye(3))) + v = [2,3,4] + rv = v.' + + @test (rv*qrmat')::RowVector == [2 3 4] +end + +@testset "Right Division" begin + mat = diagm([1,2,3]) + v = [2,3,4] + rv = v.' + + @test (rv/mat)::RowVector ≈ [2/1 3/2 4/3] + + @test (v.'/mat)::RowVector ≈ [2/1 3/2 4/3] + @test (v.'/mat.')::RowVector ≈ [2/1 3/2 4/3] + @test (rv/mat.')::RowVector ≈ [2/1 3/2 4/3] + + @test (v'/mat)::RowVector ≈ [2/1 3/2 4/3] + @test (v'/mat')::RowVector ≈ [2/1 3/2 4/3] + @test (rv/mat')::RowVector ≈ [2/1 3/2 4/3] +end + +@testset "Sparse ambiguity methods" begin + mat = sparse(diagm([1,2,3])) + v = [2,3,4] + rv = v.' + + @test (rv/mat)::RowVector ≈ [2/1 3/2 4/3] + + @test_throws Exception mat\rv +end + +@testset "AbstractTriangular ambiguity methods" begin + ut = UpperTriangular([1 0 0; 0 2 0; 0 0 3]) + v = [2,3,4] + rv = v.' + + @test (rv*ut)::RowVector == [2 6 12] + @test_throws Exception ut*rv + + @test (rv*ut.')::RowVector == [2 6 12] + @test (ut*rv.')::Vector == [2,6,12] + + @test (ut.'*rv.')::Vector == [2,6,12] + @test_throws Exception rv.'*ut.' + + @test_throws Exception ut.'*rv + @test_throws Exception rv.'*ut + + @test (rv*ut')::RowVector == [2 6 12] + @test (ut*rv')::Vector == [2,6,12] + + @test_throws Exception rv'*ut' + @test (ut'*rv')::Vector == [2,6,12] + + @test_throws Exception ut'*rv + @test_throws Exception rv'*ut + + @test (rv/ut)::RowVector ≈ [2/1 3/2 4/3] + @test (rv/ut.')::RowVector ≈ [2/1 3/2 4/3] + @test (rv/ut')::RowVector ≈ [2/1 3/2 4/3] + + @test_throws Exception ut\rv + +end + + +end # @testset "RowVector" diff --git a/test/offsetarray.jl b/test/offsetarray.jl index eb74752638473..4b408f25ac2dd 100644 --- a/test/offsetarray.jl +++ b/test/offsetarray.jl @@ -268,10 +268,11 @@ am = map(identity, a) # other functions v = OffsetArray(v0, (-3,)) @test v ≈ v -@test parent(v') == v0' -@test indices(v') === (1:1,-2:1) +@test indices(v') === (Base.OneTo(1),-2:1) A = OffsetArray(rand(4,4), (-3,5)) @test A ≈ A +@test indices(A') === (6:9, -2:1) +@test parent(A') == parent(A)' @test maximum(A) == maximum(parent(A)) @test minimum(A) == minimum(parent(A)) @test extrema(A) == extrema(parent(A)) @@ -315,10 +316,10 @@ I,J,N = findnz(z) @test mean(A_3_3) == median(A_3_3) == 5 @test mean(x->2x, A_3_3) == 10 @test mean(A_3_3, 1) == median(A_3_3, 1) == OffsetArray([2 5 8], (0,A_3_3.offsets[2])) -@test mean(A_3_3, 2) == median(A_3_3, 2) == OffsetArray([4,5,6]'', (A_3_3.offsets[1],0)) +@test mean(A_3_3, 2) == median(A_3_3, 2) == OffsetArray(reshape([4,5,6],(3,1)), (A_3_3.offsets[1],0)) @test var(A_3_3) == 7.5 @test std(A_3_3, 1) == OffsetArray([1 1 1], (0,A_3_3.offsets[2])) -@test std(A_3_3, 2) == OffsetArray([3,3,3]'', (A_3_3.offsets[1],0)) +@test std(A_3_3, 2) == OffsetArray(reshape([3,3,3], (3,1)), (A_3_3.offsets[1],0)) @test sum(OffsetArray(ones(Int,3000), -1000)) == 3000 @test vecnorm(v) ≈ vecnorm(parent(v)) diff --git a/test/show.jl b/test/show.jl index 4e070c9675b70..99dc2f7feef3e 100644 --- a/test/show.jl +++ b/test/show.jl @@ -420,7 +420,7 @@ end @test replstr(eye(10)) == "10×10 Array{Float64,2}:\n 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0\n 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0\n 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0\n 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0\n 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0\n 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0\n 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0" # an array too long vertically to fit on screen, and too long horizontally: @test replstr(collect(1.:100.)) == "100-element Array{Float64,1}:\n 1.0\n 2.0\n 3.0\n 4.0\n 5.0\n 6.0\n 7.0\n 8.0\n 9.0\n 10.0\n ⋮ \n 92.0\n 93.0\n 94.0\n 95.0\n 96.0\n 97.0\n 98.0\n 99.0\n 100.0" -@test replstr(collect(1.:100.)') == "1×100 Array{Float64,2}:\n 1.0 2.0 3.0 4.0 5.0 6.0 7.0 … 95.0 96.0 97.0 98.0 99.0 100.0" +@test replstr(collect(1.:100.)') == "1×100 RowVector{Float64,Array{Float64,1}}:\n 1.0 2.0 3.0 4.0 5.0 6.0 7.0 … 95.0 96.0 97.0 98.0 99.0 100.0" # too big in both directions to fit on screen: @test replstr((1.:100.)*(1:100)') == "100×100 Array{Float64,2}:\n 1.0 2.0 3.0 4.0 5.0 6.0 … 97.0 98.0 99.0 100.0\n 2.0 4.0 6.0 8.0 10.0 12.0 194.0 196.0 198.0 200.0\n 3.0 6.0 9.0 12.0 15.0 18.0 291.0 294.0 297.0 300.0\n 4.0 8.0 12.0 16.0 20.0 24.0 388.0 392.0 396.0 400.0\n 5.0 10.0 15.0 20.0 25.0 30.0 485.0 490.0 495.0 500.0\n 6.0 12.0 18.0 24.0 30.0 36.0 … 582.0 588.0 594.0 600.0\n 7.0 14.0 21.0 28.0 35.0 42.0 679.0 686.0 693.0 700.0\n 8.0 16.0 24.0 32.0 40.0 48.0 776.0 784.0 792.0 800.0\n 9.0 18.0 27.0 36.0 45.0 54.0 873.0 882.0 891.0 900.0\n 10.0 20.0 30.0 40.0 50.0 60.0 970.0 980.0 990.0 1000.0\n ⋮ ⋮ ⋱ \n 92.0 184.0 276.0 368.0 460.0 552.0 8924.0 9016.0 9108.0 9200.0\n 93.0 186.0 279.0 372.0 465.0 558.0 9021.0 9114.0 9207.0 9300.0\n 94.0 188.0 282.0 376.0 470.0 564.0 9118.0 9212.0 9306.0 9400.0\n 95.0 190.0 285.0 380.0 475.0 570.0 9215.0 9310.0 9405.0 9500.0\n 96.0 192.0 288.0 384.0 480.0 576.0 … 9312.0 9408.0 9504.0 9600.0\n 97.0 194.0 291.0 388.0 485.0 582.0 9409.0 9506.0 9603.0 9700.0\n 98.0 196.0 294.0 392.0 490.0 588.0 9506.0 9604.0 9702.0 9800.0\n 99.0 198.0 297.0 396.0 495.0 594.0 9603.0 9702.0 9801.0 9900.0\n 100.0 200.0 300.0 400.0 500.0 600.0 9700.0 9800.0 9900.0 10000.0" diff --git a/test/sparse/sparse.jl b/test/sparse/sparse.jl index 65f987f61351c..e3ffbf6ce9fd3 100644 --- a/test/sparse/sparse.jl +++ b/test/sparse/sparse.jl @@ -494,7 +494,7 @@ end end @testset "issue #5853, sparse diff" begin - for i=1:2, a=Any[[1 2 3], [1 2 3]', eye(3)] + for i=1:2, a=Any[[1 2 3], reshape([1, 2, 3],(3,1)), eye(3)] @test all(diff(sparse(a),i) == diff(a,i)) end end @@ -1641,9 +1641,10 @@ end end # Test temporary fix for issue #16548 in PR #16979. Brittle. Expect to remove with `\` revisions. -@testset "issue #16548" begin - @test which(\, (SparseMatrixCSC, AbstractVecOrMat)).module == Base.SparseArrays -end +# This is broken by the introduction of RowVector... see brittle comment above. +#@testset "issue #16548" begin +# @test which(\, (SparseMatrixCSC, AbstractVecOrMat)).module == Base.SparseArrays +#end @testset "row indexing a SparseMatrixCSC with non-Int integer type" begin A = sparse(UInt32[1,2,3], UInt32[1,2,3], [1.0,2.0,3.0])