From f056e29b99935cc8fbec2271c02d03e27d13c08f Mon Sep 17 00:00:00 2001 From: Klaus Crusius Date: Tue, 8 Jan 2019 16:31:58 +0100 Subject: [PATCH] proper diagonal in copytri! (fix #30055) (#30066) * proper diagonal in copytri! (fix #30055) * added sprandn methods with Type * additional parameter in copytri! for diagonal * @inline copytri! to enforce constant propagation (cherry picked from commit 4be93390b26780696cd32f262ad043ee6fbdce70) --- stdlib/LinearAlgebra/src/matmul.jl | 8 +++++--- stdlib/LinearAlgebra/src/symmetric.jl | 4 ++-- stdlib/LinearAlgebra/src/triangular.jl | 16 ++++++++-------- stdlib/LinearAlgebra/test/matmul.jl | 21 +++++++++++++++++++++ 4 files changed, 36 insertions(+), 13 deletions(-) diff --git a/stdlib/LinearAlgebra/src/matmul.jl b/stdlib/LinearAlgebra/src/matmul.jl index 9cc72742e53ff..91d72295b5422 100644 --- a/stdlib/LinearAlgebra/src/matmul.jl +++ b/stdlib/LinearAlgebra/src/matmul.jl @@ -328,14 +328,16 @@ function mul!(C::AbstractMatrix, adjA::Adjoint{<:Any,<:AbstractVecOrMat}, transB end # Supporting functions for matrix multiplication -function copytri!(A::AbstractMatrix, uplo::AbstractChar, conjugate::Bool=false) +# copy transposed(adjoint) of upper(lower) side-digonals. Optionally include diagonal. +@inline function copytri!(A::AbstractMatrix, uplo::AbstractChar, conjugate::Bool=false, diag::Bool=false) n = checksquare(A) + off = diag ? 0 : 1 if uplo == 'U' - for i = 1:(n-1), j = (i+1):n + for i = 1:n, j = (i+off):n A[j,i] = conjugate ? adjoint(A[i,j]) : transpose(A[i,j]) end elseif uplo == 'L' - for i = 1:(n-1), j = (i+1):n + for i = 1:n, j = (i+off):n A[i,j] = conjugate ? adjoint(A[j,i]) : transpose(A[j,i]) end else diff --git a/stdlib/LinearAlgebra/src/symmetric.jl b/stdlib/LinearAlgebra/src/symmetric.jl index a43e338670ea0..627a09c2a959e 100644 --- a/stdlib/LinearAlgebra/src/symmetric.jl +++ b/stdlib/LinearAlgebra/src/symmetric.jl @@ -236,14 +236,14 @@ similar(A::Union{Symmetric,Hermitian}, ::Type{T}, dims::Dims{N}) where {T,N} = s function Matrix(A::Symmetric) B = copytri!(convert(Matrix, copy(A.data)), A.uplo) for i = 1:size(A, 1) - B[i,i] = symmetric(B[i,i], sym_uplo(A.uplo))::symmetric_type(eltype(A.data)) + B[i,i] = symmetric(A[i,i], sym_uplo(A.uplo))::symmetric_type(eltype(A.data)) end return B end function Matrix(A::Hermitian) B = copytri!(convert(Matrix, copy(A.data)), A.uplo, true) for i = 1:size(A, 1) - B[i,i] = hermitian(B[i,i], sym_uplo(A.uplo))::hermitian_type(eltype(A.data)) + B[i,i] = hermitian(A[i,i], sym_uplo(A.uplo))::hermitian_type(eltype(A.data)) end return B end diff --git a/stdlib/LinearAlgebra/src/triangular.jl b/stdlib/LinearAlgebra/src/triangular.jl index 93f54f7ee5a2f..fdf135a4e0fa8 100644 --- a/stdlib/LinearAlgebra/src/triangular.jl +++ b/stdlib/LinearAlgebra/src/triangular.jl @@ -346,14 +346,14 @@ Base.copy(A::Transpose{<:Any,<:UpperTriangular}) = transpose!(copy(A.parent)) Base.copy(A::Transpose{<:Any,<:UnitLowerTriangular}) = transpose!(copy(A.parent)) Base.copy(A::Transpose{<:Any,<:UnitUpperTriangular}) = transpose!(copy(A.parent)) -transpose!(A::LowerTriangular) = UpperTriangular(copytri!(A.data, 'L')) -transpose!(A::UnitLowerTriangular) = UnitUpperTriangular(copytri!(A.data, 'L')) -transpose!(A::UpperTriangular) = LowerTriangular(copytri!(A.data, 'U')) -transpose!(A::UnitUpperTriangular) = UnitLowerTriangular(copytri!(A.data, 'U')) -adjoint!(A::LowerTriangular) = UpperTriangular(copytri!(A.data, 'L' , true)) -adjoint!(A::UnitLowerTriangular) = UnitUpperTriangular(copytri!(A.data, 'L' , true)) -adjoint!(A::UpperTriangular) = LowerTriangular(copytri!(A.data, 'U' , true)) -adjoint!(A::UnitUpperTriangular) = UnitLowerTriangular(copytri!(A.data, 'U' , true)) +transpose!(A::LowerTriangular) = UpperTriangular(copytri!(A.data, 'L', false, true)) +transpose!(A::UnitLowerTriangular) = UnitUpperTriangular(copytri!(A.data, 'L', false, true)) +transpose!(A::UpperTriangular) = LowerTriangular(copytri!(A.data, 'U', false, true)) +transpose!(A::UnitUpperTriangular) = UnitLowerTriangular(copytri!(A.data, 'U', false, true)) +adjoint!(A::LowerTriangular) = UpperTriangular(copytri!(A.data, 'L' , true, true)) +adjoint!(A::UnitLowerTriangular) = UnitUpperTriangular(copytri!(A.data, 'L' , true, true)) +adjoint!(A::UpperTriangular) = LowerTriangular(copytri!(A.data, 'U' , true, true)) +adjoint!(A::UnitUpperTriangular) = UnitLowerTriangular(copytri!(A.data, 'U' , true, true)) diag(A::LowerTriangular) = diag(A.data) diag(A::UnitLowerTriangular) = fill(one(eltype(A)), size(A,1)) diff --git a/stdlib/LinearAlgebra/test/matmul.jl b/stdlib/LinearAlgebra/test/matmul.jl index 845a0b5598169..bc3d3d8607afa 100644 --- a/stdlib/LinearAlgebra/test/matmul.jl +++ b/stdlib/LinearAlgebra/test/matmul.jl @@ -301,6 +301,27 @@ end @test_throws ArgumentError LinearAlgebra.copytri!(Matrix{Float64}(undef,10,10),'Z') +@testset "Issue 30055" begin + B = [1+im 2+im 3+im; 4+im 5+im 6+im; 7+im 9+im im] + A = UpperTriangular(B) + @test copy(transpose(A)) == transpose(A) + @test copy(A') == A' + A = LowerTriangular(B) + @test copy(transpose(A)) == transpose(A) + @test copy(A') == A' + B = Matrix{Matrix{Complex{Int}}}(undef, 2, 2) + B[1,1] = [1+im 2+im; 3+im 4+im] + B[2,1] = [1+2im 1+3im;1+3im 1+4im] + B[1,2] = [7+im 8+2im; 9+3im 4im] + B[2,2] = [9+im 8+im; 7+im 6+im] + A = UpperTriangular(B) + @test copy(transpose(A)) == transpose(A) + @test copy(A') == A' + A = LowerTriangular(B) + @test copy(transpose(A)) == transpose(A) + @test copy(A') == A' +end + @testset "gemv! and gemm_wrapper for $elty" for elty in [Float32,Float64,ComplexF64,ComplexF32] A10x10, x10, x11 = Array{elty}.(undef, ((10,10), 10, 11)) @test_throws DimensionMismatch LinearAlgebra.gemv!(x10,'N',A10x10,x11)