Skip to content

Commit

Permalink
Work on implementing non-recursive transpose.
Browse files Browse the repository at this point in the history
 * RowVector no longer transposes elements
 * MappedArray used to propagate conj, transpose, adjoint, etc
  • Loading branch information
Andy Ferris committed Nov 22, 2017
1 parent ea34968 commit 6c60c1e
Show file tree
Hide file tree
Showing 26 changed files with 532 additions and 400 deletions.
8 changes: 8 additions & 0 deletions NEWS.md
Original file line number Diff line number Diff line change
Expand Up @@ -290,6 +290,14 @@ This section lists changes that do not have deprecation warnings.
Its return value has been removed. Use the `process_running` function
to determine if a process has already exited.

* `transpose` and `transpose!` no longer recursively transpose the elements of the
container. Similarly, `RowVector` no longer provides a transposed view of the elements.
Transposition now simply rearranges the elements of containers of data, such as arrays
of strings. Note that the renamed `adjoint` method (formerly `ctranspose`) does still
act in a recursive manner, and that (very occassionally) `conj(adjoint(...))` will be
preferrable to `transpose` for linear algebra problems using nested arrays as "block
matrices". ([#23424])

Library improvements
--------------------

Expand Down
88 changes: 88 additions & 0 deletions base/abstractarraymath.jl
Original file line number Diff line number Diff line change
Expand Up @@ -83,6 +83,94 @@ end

squeeze(A::AbstractArray, dim::Integer) = squeeze(A, (Int(dim),))

## Transposition ##

"""
transpose(A::AbstractMatrix)
The transposition operator (`.'`).
# Examples
```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

transpose(a::AbstractArray) = error("transpose not defined for $(typeof(a)). Consider using `permutedims` for higher-dimensional arrays.")

"""
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!(identity, 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

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

## Unary operators ##

Expand Down
3 changes: 3 additions & 0 deletions base/exports.jl
Original file line number Diff line number Diff line change
Expand Up @@ -74,6 +74,9 @@ export
LinSpace,
LowerTriangular,
Irrational,
MappedArray,
MappedVector,
MappedMatrix,
Matrix,
MergeSort,
NTuple,
Expand Down
124 changes: 124 additions & 0 deletions base/linalg/adjoint.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,124 @@
# This file is a part of Julia. License is MIT: https://julialang.org/license

adjoint(a::AbstractArray) = error("adjoint not defined for $(typeof(a)). Consider using `permutedims` for higher-dimensional arrays.")

## adjoint ##

"""
adjoint(v::AbstractVector)
Creates a `RowVector` from `v` where `adjoint` has been applied recursively to the elements.
Conceptually, this is intended create the "dual vector" of `v` (note however that the output
is strictly an `AbstractMatrix`). Note also that the output is a view of `v`.
"""
@inline adjoint(vec::AbstractVector) = RowVector(_map(adjoint, vec))
@inline adjoint(rowvec::RowVector) = _map(adjoint, parent(rowvec))

"""
adjoint(m::AbstractMatrix)
Returns the Hermitian adjoint of `m`, where `m` has been transposed and `adjoint` is applied
recursively to the elements.
"""
function adjoint(a::AbstractMatrix)
(ind1, ind2) = indices(a)
b = similar(a, promote_op(adjoint, eltype(a)), (ind2, ind1))
adjoint!(b, a)
end

"""
adjoint!(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.
"""
adjoint!(b::AbstractMatrix, a::AbstractMatrix) = transpose_f!(adjoint, b, a)
function adjoint!(b::AbstractVector, a::AbstractMatrix)
if indices(b, 1) != indices(a, 2) || indices(a, 1) != 1:1
throw(DimensionMismatch("adjoint"))
end
adjointcopy!(b, a)
end
function adjoint!(b::AbstractMatrix, a::AbstractVector)
if indices(b, 2) != indices(a, 1) || indices(b, 1) != 1:1
throw(DimensionMismatch("adjoint"))
end
adjointcopy!(b, a)
end

function adjointcopy!(b, a)
ra = eachindex(a)
rb = eachindex(b)
if rb == ra
for i rb
b[i] = adjoint(a[i])
end
else
for (i, j) zip(rb, ra)
b[i] = adjoint(a[j])
end
end
end

"""
adjoint(x::Number)
The (complex) conjugate of `x`, `conj(x)`.
"""
adjoint(x::Number) = conj(x)


## conjadjoint ##

"""
conjadjoint(a)
Returns `conj(adjoint(a))`.
"""
conjadjoint(a::AbstractArray) = conj(adjoint(a))

(::typeof(conj), ::typeof(conj)) = identity
(::typeof(adjoint), ::typeof(adjoint)) = identity
(::typeof(conjadjoint), ::typeof(conjadjoint)) = identity

(::typeof(conj), ::typeof(adjoint)) = conjadjoint
(::typeof(adjoint), ::typeof(conj)) = conjadjoint
(::typeof(conj), ::typeof(conjadjoint)) = adjoint
(::typeof(conjadjoint), ::typeof(conj)) = adjoint
(::typeof(adjoint), ::typeof(conjadjoint)) = conj
(::typeof(conjadjoint), ::typeof(adjoint)) = conj

## mapped array aliases ##

const ConjArray{T,N,A<:AbstractArray{<:Any,N}} = MappedArray{T,N,typeof(conj),typeof(conj),A}
const ConjVector{T,A<:AbstractVector} = MappedVector{T,typeof(conj),typeof(conj),A}
const ConjMatrix{T,A<:AbstractMatrix} = MappedMatrix{T,typeof(conj),typeof(conj),A}

inv_func(::typeof(adjoint)) = adjoint
inv_func(::typeof(conjadjoint)) = conjadjoint

@inline _map(f, a::AbstractArray) = MappedArray(f, a)
@inline _map(f, a::MappedArray) = map(f, a)

## lazy conj ##

"""
conj(v::RowVector)
Returns a [`ConjArray`](@ref) lazy view of the input, where each element is conjugated.
# Examples
```jldoctest
julia> v = [1+im, 1-im].'
1×2 RowVector{Complex{Int64},Array{Complex{Int64},1}}:
1+1im 1-1im
julia> conj(v)
1×2 RowVector{Complex{Int64},ConjArray{Complex{Int64},1,Array{Complex{Int64},1}}}:
1-1im 1+1im
```
"""
@inline conj(rowvec::RowVector) = RowVector(_map(conj, parent(rowvec)))

3 changes: 2 additions & 1 deletion base/linalg/bidiag.jl
Original file line number Diff line number Diff line change
Expand Up @@ -243,7 +243,8 @@ broadcast(::typeof(floor), ::Type{T}, M::Bidiagonal) where {T<:Integer} = Bidiag
broadcast(::typeof(ceil), ::Type{T}, M::Bidiagonal) where {T<:Integer} = Bidiagonal(ceil.(T, M.dv), ceil.(T, M.ev), M.uplo)

transpose(M::Bidiagonal) = Bidiagonal(M.dv, M.ev, M.uplo == 'U' ? :L : :U)
adjoint(M::Bidiagonal) = Bidiagonal(conj(M.dv), conj(M.ev), M.uplo == 'U' ? :L : :U)
adjoint(M::Bidiagonal{<:Number}) = Bidiagonal(conj(M.dv), conj(M.ev), M.uplo == 'U' ? :L : :U)
adjoint(M::Bidiagonal) = Bidiagonal(adjoint.(M.dv), adjoint.(M.ev), M.uplo == 'U' ? :L : :U)

istriu(M::Bidiagonal) = M.uplo == 'U' || iszero(M.ev)
istril(M::Bidiagonal) = M.uplo == 'L' || iszero(M.ev)
Expand Down
61 changes: 0 additions & 61 deletions base/linalg/conjarray.jl

This file was deleted.

15 changes: 7 additions & 8 deletions base/linalg/diagonal.jl
Original file line number Diff line number Diff line change
Expand Up @@ -197,7 +197,7 @@ function Ac_mul_B(A::AbstractMatrix, D::Diagonal)
A_mul_B!(Ac, D)
end

At_mul_B(D::Diagonal, B::Diagonal) = Diagonal(transpose.(D.diag) .* B.diag)
At_mul_B(D::Diagonal, B::Diagonal) = Diagonal(D.diag .* B.diag)
At_mul_B(A::AbstractTriangular, D::Diagonal) = A_mul_B!(transpose(A), D)
function At_mul_B(A::AbstractMatrix, D::Diagonal)
At = similar(A, promote_op(*, eltype(A), eltype(D.diag)), (size(A, 2), size(A, 1)))
Expand All @@ -214,7 +214,7 @@ function A_mul_Bc(D::Diagonal, A::AbstractMatrix)
A_mul_B!(D, Ac)
end

A_mul_Bt(D::Diagonal, B::Diagonal) = Diagonal(D.diag .* transpose.(B.diag))
A_mul_Bt(D::Diagonal, B::Diagonal) = Diagonal(D.diag .* B.diag)
A_mul_Bt(D::Diagonal, B::AbstractTriangular) = A_mul_B!(D, transpose(B))
function A_mul_Bt(D::Diagonal, A::AbstractMatrix)
At = similar(A, promote_op(*, eltype(A), eltype(D.diag)), (size(A, 2), size(A, 1)))
Expand All @@ -223,7 +223,7 @@ function A_mul_Bt(D::Diagonal, A::AbstractMatrix)
end

Ac_mul_Bc(D::Diagonal, B::Diagonal) = Diagonal(adjoint.(D.diag) .* adjoint.(B.diag))
At_mul_Bt(D::Diagonal, B::Diagonal) = Diagonal(transpose.(D.diag) .* transpose.(B.diag))
At_mul_Bt(D::Diagonal, B::Diagonal) = Diagonal(D.diag .* B.diag)

A_mul_B!(A::Diagonal,B::Diagonal) = throw(MethodError(A_mul_B!, Tuple{Diagonal,Diagonal}))
At_mul_B!(A::Diagonal,B::Diagonal) = throw(MethodError(At_mul_B!, Tuple{Diagonal,Diagonal}))
Expand All @@ -239,14 +239,14 @@ A_mul_Bc!(A::AbstractMatrix,B::Diagonal) = scale!(A,conj(B.diag))
# Get ambiguous method if try to unify AbstractVector/AbstractMatrix here using AbstractVecOrMat
A_mul_B!(out::AbstractVector, A::Diagonal, in::AbstractVector) = out .= A.diag .* in
Ac_mul_B!(out::AbstractVector, A::Diagonal, in::AbstractVector) = out .= adjoint.(A.diag) .* in
At_mul_B!(out::AbstractVector, A::Diagonal, in::AbstractVector) = out .= transpose.(A.diag) .* in
At_mul_B!(out::AbstractVector, A::Diagonal, in::AbstractVector) = out .= A.diag .* in

A_mul_B!(out::AbstractMatrix, A::Diagonal, in::AbstractMatrix) = out .= A.diag .* in
Ac_mul_B!(out::AbstractMatrix, A::Diagonal, in::AbstractMatrix) = out .= adjoint.(A.diag) .* in
At_mul_B!(out::AbstractMatrix, A::Diagonal, in::AbstractMatrix) = out .= transpose.(A.diag) .* in
At_mul_B!(out::AbstractMatrix, A::Diagonal, in::AbstractMatrix) = out .= A.diag .* in

# ambiguities with Symmetric/Hermitian
# RealHermSymComplex[Sym]/[Herm] only include Number; invariant to [c]transpose
# RealHermSymComplex[Sym]/[Herm] only include Number; invariant to transpose/adjoint
A_mul_Bt(A::Diagonal, B::RealHermSymComplexSym) = A*B
At_mul_B(A::RealHermSymComplexSym, B::Diagonal) = A*B
A_mul_Bc(A::Diagonal, B::RealHermSymComplexHerm) = A*B
Expand Down Expand Up @@ -317,8 +317,7 @@ Ac_ldiv_B(F::Factorization, D::Diagonal) =
@inline A_mul_Bc(D::Diagonal, rowvec::RowVector) = D*adjoint(rowvec)

conj(D::Diagonal) = Diagonal(conj(D.diag))
transpose(D::Diagonal{<:Number}) = D
transpose(D::Diagonal) = Diagonal(transpose.(D.diag))
transpose(D::Diagonal) = D
adjoint(D::Diagonal{<:Number}) = conj(D)
adjoint(D::Diagonal) = Diagonal(adjoint.(D.diag))

Expand Down
2 changes: 1 addition & 1 deletion base/linalg/factorization.jl
Original file line number Diff line number Diff line change
Expand Up @@ -81,7 +81,7 @@ end

# fallback methods for transposed solves
At_ldiv_B(F::Factorization{<:Real}, B::AbstractVecOrMat) = Ac_ldiv_B(F, B)
At_ldiv_B(F::Factorization, B) = conj.(Ac_ldiv_B(F, conj.(B)))
At_ldiv_B(F::Factorization, B) = conj.(Ac_ldiv_B(F, conj.(B))) # TODO fix for array eltypes?

"""
A_ldiv_B!([Y,] A, B) -> Y
Expand Down
Loading

0 comments on commit 6c60c1e

Please sign in to comment.