Skip to content

Commit

Permalink
Introduce RowVector as the transpose of a vector (#19670)
Browse files Browse the repository at this point in the history 
* Moved all array transpose functions to LinAlg

Transposition is a concept of linear algebra rather than
multidimensional arrays of data.

* 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
  • Loading branch information
@andyferris @stevengj
andyferris authored and stevengj committed Jan 13, 2017
1 parent aa2ea55 commit 3d6fa35
Show file tree
Hide file tree
Showing 32 changed files with 833 additions and 283 deletions.
4 changes: 4 additions & 0 deletions NEWS.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
--------------------

Expand Down
23 changes: 0 additions & 23 deletions base/abstractarray.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -644,29 +644,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)
Expand Down
2 changes: 0 additions & 2 deletions base/abstractarraymath.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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 ##

Expand Down
127 changes: 0 additions & 127 deletions base/arraymath.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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 ]
88 changes: 0 additions & 88 deletions base/bitarray.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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 ##

Expand Down
10 changes: 10 additions & 0 deletions base/docs/basedocs.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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".'"

Expand Down
1 change: 1 addition & 0 deletions base/exports.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -98,6 +98,7 @@ export
RoundNearestTiesUp,
RoundToZero,
RoundUp,
RowVector,
AbstractSerializer,
SerializationState,
Set,
Expand Down
9 changes: 9 additions & 0 deletions base/linalg/bidiag.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
Loading

0 comments on commit 3d6fa35

Please sign in to comment.