From db90a8e5e89189fce33b3c5d63c769c75b18b2aa Mon Sep 17 00:00:00 2001 From: Andy Ferris Date: Sun, 15 Jan 2017 03:22:16 -0800 Subject: [PATCH] Added `ConjArray` wrapper type for conjugate views By default, this is used only in conjugation with `RowVector`, so that both `transpose(vec)` and `ctranspose(vec)` both return views. --- NEWS.md | 4 +++ base/exports.jl | 3 ++ base/linalg/conjarray.jl | 62 ++++++++++++++++++++++++++++++++++++++++ base/linalg/linalg.jl | 4 +++ base/linalg/rowvector.jl | 62 +++++++++++++++++++++++++++------------- doc/src/stdlib/linalg.md | 1 + test/choosetests.jl | 2 +- test/linalg/conjarray.jl | 23 +++++++++++++++ test/linalg/rowvector.jl | 10 ++----- test/sparse/sparse.jl | 10 +++---- 10 files changed, 148 insertions(+), 33 deletions(-) create mode 100644 base/linalg/conjarray.jl create mode 100644 test/linalg/conjarray.jl diff --git a/NEWS.md b/NEWS.md index a7dc330a5e45e..2d752f7bed7b3 100644 --- a/NEWS.md +++ b/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 ----------------------------- diff --git a/base/exports.jl b/base/exports.jl index 78b9ee477cc2b..c044051ddb263 100644 --- a/base/exports.jl +++ b/base/exports.jl @@ -50,6 +50,9 @@ export Complex128, Complex64, Complex32, + ConjArray, + ConjVector, + ConjMatrix, DenseMatrix, DenseVecOrMat, DenseVector, diff --git a/base/linalg/conjarray.jl b/base/linalg/conjarray.jl new file mode 100644 index 0000000000000..d836994cfda84 --- /dev/null +++ b/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) diff --git a/base/linalg/linalg.jl b/base/linalg/linalg.jl index 9ccc8ae4c827a..e42f68cae28be 100644 --- a/base/linalg/linalg.jl +++ b/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") diff --git a/base/linalg/rowvector.jl b/base/linalg/rowvector.jl index 980363d3ea5db..b8d49eabc077b 100644 --- a/base/linalg/rowvector.jl +++ b/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))")) diff --git a/doc/src/stdlib/linalg.md b/doc/src/stdlib/linalg.md index a3a3a74cb8bd9..f8dc0fe81446d 100644 --- a/doc/src/stdlib/linalg.md +++ b/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 diff --git a/test/choosetests.jl b/test/choosetests.jl index 67f32c8d0d190..a81b1f9d98ff6 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/rowvector"] + "linalg/hessenberg", "linalg/rowvector", "linalg/conjarray"] if Base.USE_GPL_LIBS push!(linalgtests, "linalg/arnoldi") end diff --git a/test/linalg/conjarray.jl b/test/linalg/conjarray.jl new file mode 100644 index 0000000000000..7f3f4190b2ad9 --- /dev/null +++ b/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 diff --git a/test/linalg/rowvector.jl b/test/linalg/rowvector.jl index e2efb7c189c95..f188a3082c000 100644 --- a/test/linalg/rowvector.jl +++ b/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,7 +102,7 @@ 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 @@ -112,7 +110,7 @@ end @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" diff --git a/test/sparse/sparse.jl b/test/sparse/sparse.jl index f6bed5fe0dfe3..baec169e41449 100644 --- a/test/sparse/sparse.jl +++ b/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])