Added ConjArray wrapper type for conjugate array views

NEWS.md

@@ -282,6 +282,10 @@ Library improvements
 
   * New `@macroexpand` macro as a convenient alternative to the `macroexpand` function ([#18660]).
 
+  * Introduced a wrapper type for lazy complex conjugation of arrays, `ConjArray`.
+    Currently, it is used by default for the new `RowVector` type only, and
+    enforces that both `transpose(vec)` and `ctranspose(vec)` are views not copies ([#20047]).
+
 Compiler/Runtime improvements
 -----------------------------
 

base/exports.jl

@@ -50,6 +50,9 @@ export
     Complex128,
     Complex64,
     Complex32,
+    ConjArray,
+    ConjVector,
+    ConjMatrix,
     DenseMatrix,
     DenseVecOrMat,
     DenseVector,

base/linalg/conjarray.jl

@@ -0,0 +1,62 @@
+# This file is a part of Julia. License is MIT: http://julialang.org/license
+
+"""
+    ConjArray(array)
+
+A lazy-view wrapper of an `AbstractArray`, taking the elementwise complex conjugate. This
+type is usually constructed (and unwrapped) via the [`conj`](@ref) function (or related
+[`ctranspose`](@ref)), but currently this is the default behavior for `RowVector` only. For
+other arrays, the `ConjArray` constructor can be used directly.
+
+# Examples
+
+```jldoctest
+julia> [1+im, 1-im]'
+1×2 RowVector{Complex{Int64},ConjArray{Complex{Int64},1,Array{Complex{Int64},1}}}:
+ 1-1im  1+1im
+
+julia> ConjArray([1+im 0; 0 1-im])
+2×2 ConjArray{Complex{Int64},2,Array{Complex{Int64},2}}:
+ 1-1im  0+0im
+ 0+0im  1+1im
+```
+"""
+immutable ConjArray{T, N, A <: AbstractArray} <: AbstractArray{T, N}
+    parent::A
+end
+
+@inline ConjArray{T,N}(a::AbstractArray{T,N}) = ConjArray{conj_type(T), N, typeof(a)}(a)
+
+typealias ConjVector{T, V <: AbstractVector} ConjArray{T, 1, V}
+@inline ConjVector{T}(v::AbstractVector{T}) = ConjArray{conj_type(T), 1, typeof(v)}(v)
+
+typealias ConjMatrix{T, M <: AbstractMatrix} ConjArray{T, 2, M}
+@inline ConjMatrix{T}(m::AbstractMatrix{T}) = ConjArray{conj_type(T), 2, typeof(m)}(m)
+
+# This type can cause the element type to change under conjugation - e.g. an array of complex arrays.
+@inline conj_type(x) = conj_type(typeof(x))
+@inline conj_type{T}(::Type{T}) = promote_op(conj, T)
+
+@inline parent(c::ConjArray) = c.parent
+@inline parent_type(c::ConjArray) = parent_type(typeof(c))
+@inline parent_type{T,N,A}(::Type{ConjArray{T,N,A}}) = A
+
+@inline size(a::ConjArray) = size(a.parent)
+linearindexing{CA <: ConjArray}(::CA) = linearindexing(parent_type(CA))
+linearindexing{CA <: ConjArray}(::Type{CA}) = linearindexing(parent_type(CA))
+
+@propagate_inbounds getindex{T,N}(a::ConjArray{T,N}, i::Int) = conj(getindex(a.parent, i))
+@propagate_inbounds getindex{T,N}(a::ConjArray{T,N}, i::Vararg{Int,N}) = conj(getindex(a.parent, i...))
+@propagate_inbounds setindex!{T,N}(a::ConjArray{T,N}, v, i::Int) = setindex!(a.parent, conj(v), i)
+@propagate_inbounds setindex!{T,N}(a::ConjArray{T,N}, v, i::Vararg{Int,N}) = setindex!(a.parent, conj(v), i...)
+
+@inline similar{T,N}(a::ConjArray, ::Type{T}, dims::Dims{N}) = similar(parent(a), T, dims)
+
+# Currently, this is default behavior for RowVector only
+@inline conj(a::ConjArray) = parent(a)
+
+# Helper functions, currently used by RowVector
+@inline _conj(a::AbstractArray) = ConjArray(a)
+@inline _conj{T<:Real}(a::AbstractArray{T}) = a
+@inline _conj(a::ConjArray) = parent(a)
+@inline _conj{T<:Real}(a::ConjArray{T}) = parent(a)

base/linalg/linalg.jl

@@ -22,6 +22,9 @@ export
 
 # Types
     RowVector,
+    ConjArray,
+    ConjVector,
+    ConjMatrix,
     SymTridiagonal,
     Tridiagonal,
     Bidiagonal,
@@ -237,6 +240,7 @@ 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("conjarray.jl")
 include("transpose.jl")
 include("rowvector.jl")
 

base/linalg/rowvector.jl

@@ -1,17 +1,15 @@
 """
     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.
+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`](@ref)
+function or `.'` operator (or related [`ctranspose`](@ref) 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
@@ -21,11 +19,12 @@ immutable RowVector{T,V<:AbstractVector} <: AbstractMatrix{T}
     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}`")
 
+typealias ConjRowVector{T, CV <: ConjVector} RowVector{T, CV}
+
 # The element type may be transformed as transpose is recursive
 @inline transpose_type{T}(::Type{T}) = promote_op(transpose, T)
 
@@ -47,10 +46,12 @@ end
 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
+# similar tries to maintain the RowVector wrapper and the parent type
+@inline similar(rowvec::RowVector) = RowVector(similar(parent(rowvec)))
+@inline similar{T}(rowvec::RowVector, ::Type{T}) = RowVector(similar(parent(rowvec), transpose_type(T)))
+
+# Resizing similar currently loses its RowVector property.
+@inline similar{T,N}(rowvec::RowVector, ::Type{T}, dims::Dims{N}) = similar(parent(rowvec), T, dims)
 
 # Basic methods
 """
@@ -73,18 +74,34 @@ julia> transpose(v)
 ```
 """
 @inline transpose(vec::AbstractVector) = RowVector(vec)
-@inline ctranspose{T}(vec::AbstractVector{T}) = RowVector(conj(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 transpose(rowvec::ConjRowVector) = copy(rowvec.vec) # remove the ConjArray wrapper from any raw vector
 @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))
+"""
+    conj(rowvector)
+
+Returns a [`ConjArray`](@ref) lazy view of the input, where each element is conjugated.
+
+### Example
+
+```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(_conj(rowvec.vec))
 @inline conj{T<:Real}(rowvec::RowVector{T}) = rowvec
 
 # AbstractArray interface
@@ -147,6 +164,11 @@ end
 
 # Multiplication #
 
+# inner product -> dot product specializations
+@inline *{T<:Real}(rowvec::RowVector{T}, vec::AbstractVector{T}) = dot(parent(rowvec), vec)
+@inline *(rowvec::ConjRowVector, vec::AbstractVector) = dot(rowvec', vec)
+
+# Generic behavior
 @inline function *(rowvec::RowVector, vec::AbstractVector)
     if length(rowvec) != length(vec)
         throw(DimensionMismatch("A has dimensions $(size(rowvec)) but B has dimensions $(size(vec))"))

doc/src/stdlib/linalg.md

@@ -100,6 +100,7 @@ Base.LinAlg.istriu
 Base.LinAlg.isdiag
 Base.LinAlg.ishermitian
 Base.LinAlg.RowVector
+Base.LinAlg.ConjArray
 Base.transpose
 Base.transpose!
 Base.ctranspose

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/rowvector"]
+                   "linalg/hessenberg", "linalg/rowvector", "linalg/conjarray"]
     if Base.USE_GPL_LIBS
         push!(linalgtests, "linalg/arnoldi")
     end

test/linalg/conjarray.jl

@@ -0,0 +1,23 @@
+# This file is a part of Julia. License is MIT: http://julialang.org/license
+
+@testset "Core" begin
+    m = [1+im 2; 2 4-im]
+    cm = ConjArray(m)
+    @test cm[1,1] == 1-im
+    @test trace(cm*m) == 27
+
+    v = [[1+im], [1-im]]
+    cv = ConjArray(v)
+    @test cv[1] == [1-im]
+end
+
+@testset "RowVector conjugates" begin
+    v = [1+im, 1-im]
+    rv = v'
+    @test (parent(rv) isa ConjArray)
+    @test rv' === v
+
+    # Currently, view behavior defaults to only RowVectors.
+    @test isa((v').', Vector)
+    @test isa((v.')', Vector)
+end

test/linalg/rowvector.jl

@@ -1,7 +1,5 @@
 # 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]
@@ -104,15 +102,15 @@ end
 
     @test (rv*v) === 14
     @test (rv*mat)::RowVector == [1 4 9]
-    @test [1]*reshape([1],(1,1)) == reshape([1],(1,1))
+    @test [1]*reshape([1],(1,1)) == reshape([1], (1,1))
     @test_throws DimensionMismatch rv*rv
     @test (v*rv)::Matrix == [1 2 3; 2 4 6; 3 6 9]
     @test_throws DimensionMismatch v*v # Was previously a missing method error, now an error message
     @test_throws DimensionMismatch mat*rv
 
     @test_throws DimensionMismatch rv*v.'
     @test (rv*mat.')::RowVector == [1 4 9]
-    @test [1]*reshape([1],(1,1)).' == reshape([1],(1,1))
+    @test [1]*reshape([1],(1,1)).' == reshape([1], (1,1))
     @test rv*rv.' === 14
     @test_throws DimensionMismatch v*rv.'
     @test (v*v.')::Matrix == [1 2 3; 2 4 6; 3 6 9]
@@ -142,7 +140,7 @@ end
 
     @test_throws DimensionMismatch cz*z'
     @test (cz*mat')::RowVector == [-2im 4 9]
-    @test [1]*reshape([1],(1,1))' == reshape([1],(1,1))
+    @test [1]*reshape([1],(1,1))' == reshape([1], (1,1))
     @test cz*cz' === 15 + 0im
     @test_throws DimensionMismatch z*cz'
     @test (z*z')::Matrix == [2 2+2im 3+3im; 2-2im 4 6; 3-3im 6 9]
@@ -251,5 +249,3 @@ end
         @test A'*x' == A'*y == B*x' == B*y == C'
     end
 end
-
-end # @testset "RowVector"

test/sparse/sparse.jl

@@ -1641,11 +1641,11 @@ end
     @test At_ldiv_B(ltintmat, sparse(intmat)) ≈ At_ldiv_B(ltintmat, intmat)
 end
 
-# Test temporary fix for issue #16548 in PR #16979. Brittle. Expect to remove with `\` revisions.
-# This is broken by the introduction of RowVector... see brittle comment above.
-#@testset "issue #16548" begin
-#    @test which(\, (SparseMatrixCSC, AbstractVecOrMat)).module == Base.SparseArrays
-#end
+# Test temporary fix for issue #16548 in PR #16979. Somewhat brittle. Expect to remove with `\` revisions.
+@testset "issue #16548" begin
+    ms = methods(\, (SparseMatrixCSC, AbstractVecOrMat)).ms
+    @test all(m -> m.module == Base.SparseArrays, ms)
+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])