Make adjoint/transpose for Diagonal/Bidiagonal lazy (#1368)

Currently, `adjoint` and `transpose` for a block `Diagonal` and a block
`Bidiagonal` matrix allocates the diagonals. This was not intentional,
and we may use lazy wrappers instead that ensures that the result shares
memory with the original array, which is what the documentation states.
After this PR, mutating the `adjoint` will automatically mutate the
parent array.

Author: jishnub

src/adjtrans.jl

@@ -552,3 +552,20 @@ conj(A::Adjoint) = transpose(A.parent)
 function Base.replace_in_print_matrix(A::AdjOrTrans,i::Integer,j::Integer,s::AbstractString)
     Base.replace_in_print_matrix(parent(A), j, i, s)
 end
+
+# Special adjoint/transpose methods for Adjoint/Transpose that have been reshaped as a vector
+# these are used in transposing banded matrices by forwarding the operation to the bands
+"""
+    _vectranspose(A::AbstractVector)::AbstractVector
+
+Compute `vec(transpose(A))`, but avoid an allocating reshape if possible
+"""
+_vectranspose(A::AbstractVector) = vec(transpose(A))
+_vectranspose(A::Base.ReshapedArray{<:Any,1,<:TransposeAbsVec}) = transpose(parent(A))
+"""
+    _vecadjoint(A::AbstractVector)::AbstractVector
+
+Compute `vec(adjoint(A))`, but avoid an allocating reshape if possible
+"""
+_vecadjoint(A::AbstractVector) = vec(adjoint(A))
+_vecadjoint(A::Base.ReshapedArray{<:Any,1,<:AdjointAbsVec}) = adjoint(parent(A))

src/bidiag.jl

@@ -308,12 +308,8 @@ for func in (:conj, :copy, :real, :imag)
 end
 isreal(M::Bidiagonal) = isreal(M.dv) && isreal(M.ev)
 
-adjoint(B::Bidiagonal{<:Number}) = Bidiagonal(vec(adjoint(B.dv)), vec(adjoint(B.ev)), B.uplo == 'U' ? :L : :U)
-adjoint(B::Bidiagonal{<:Number, <:Base.ReshapedArray{<:Number,1,<:Adjoint}}) =
-    Bidiagonal(adjoint(parent(B.dv)), adjoint(parent(B.ev)), B.uplo == 'U' ? :L : :U)
-adjoint(B::Bidiagonal) = Bidiagonal(adjoint.(B.dv), adjoint.(B.ev), B.uplo == 'U' ? :L : :U)
-transpose(B::Bidiagonal{<:Number}) = Bidiagonal(B.dv, B.ev, B.uplo == 'U' ? :L : :U)
-transpose(B::Bidiagonal) = Bidiagonal(transpose.(B.dv), transpose.(B.ev), B.uplo == 'U' ? :L : :U)
+adjoint(B::Bidiagonal) = Bidiagonal(_vecadjoint(B.dv), _vecadjoint(B.ev), B.uplo == 'U' ? :L : :U)
+transpose(B::Bidiagonal) = Bidiagonal(_vectranspose(B.dv), _vectranspose(B.ev), B.uplo == 'U' ? :L : :U)
 permutedims(B::Bidiagonal) = Bidiagonal(B.dv, B.ev, B.uplo == 'U' ? 'L' : 'U')
 function permutedims(B::Bidiagonal, perm)
     Base.checkdims_perm(axes(B), axes(B), perm)

src/diagonal.jl

@@ -906,11 +906,8 @@ end
 end
 
 conj(D::Diagonal) = Diagonal(conj(D.diag))
-transpose(D::Diagonal{<:Number}) = D
-transpose(D::Diagonal) = Diagonal(transpose.(D.diag))
-adjoint(D::Diagonal{<:Number}) = Diagonal(vec(adjoint(D.diag)))
-adjoint(D::Diagonal{<:Number,<:Base.ReshapedArray{<:Number,1,<:Adjoint}}) = Diagonal(adjoint(parent(D.diag)))
-adjoint(D::Diagonal) = Diagonal(adjoint.(D.diag))
+transpose(D::Diagonal) = Diagonal(_vectranspose(D.diag))
+adjoint(D::Diagonal) = Diagonal(_vecadjoint(D.diag))
 permutedims(D::Diagonal) = D
 permutedims(D::Diagonal, perm) = (Base.checkdims_perm(axes(D), axes(D), perm); D)
 

test/bidiag.jl

@@ -1240,4 +1240,19 @@ end
     @test_throws BoundsError B[LinearAlgebra.BandIndex(0,size(B,1)+1)]
 end
 
+@testset "lazy adjtrans" begin
+    B = Bidiagonal(fill([1 2; 3 4], 3), fill([5 6; 7 8], 2), :U)
+    m = [2 4; 6 8]
+    for op in (transpose, adjoint)
+        C = op(B)
+        el = op(m)
+        C[1,1] = el
+        @test B[1,1] == m
+        C[2,1] = el
+        @test B[1,2] == m
+        @test (@allocated op(B)) == 0
+        @test (@allocated op(op(B))) == 0
+    end
+end
+
 end # module TestBidiagonal

test/diagonal.jl

@@ -1515,4 +1515,17 @@ end
     @test_throws BoundsError D[LinearAlgebra.BandIndex(0,size(D,1)+1)]
 end
 
+@testset "lazy adjtrans" begin
+    D = Diagonal(fill([1 2; 3 4], 3))
+    m = [2 4; 6 8]
+    for op in (transpose, adjoint)
+        C = op(D)
+        el = op(m)
+        C[1,1] = el
+        @test D[1,1] == m
+        @test (@allocated op(D)) == 0
+        @test (@allocated op(op(D))) == 0
+    end
+end
+
 end # module TestDiagonal