From a8f22c4cee2f6b1e43fc139fd15ff1348142d9ef Mon Sep 17 00:00:00 2001 From: Jishnu Bhattacharya Date: Thu, 10 Apr 2025 16:09:22 +0530 Subject: [PATCH 1/8] Matrix constructor for triangular --- src/generic.jl | 28 ++++++++++++++++++---------- src/triangular.jl | 6 ++++++ 2 files changed, 24 insertions(+), 10 deletions(-) diff --git a/src/generic.jl b/src/generic.jl index de37f081..4dea4643 100644 --- a/src/generic.jl +++ b/src/generic.jl @@ -2079,21 +2079,29 @@ function copytrito!(B::AbstractMatrix, A::AbstractMatrix, uplo::AbstractChar) A = Base.unalias(B, A) if uplo == 'U' LAPACK.lacpy_size_check(size(B), (n < m ? n : m, n)) - for j in axes(A,2), i in axes(A,1)[begin : min(j,end)] - # extract the parents for UpperTriangular matrices - Bv, Av = uppertridata(B), uppertridata(A) - @inbounds Bv[i,j] = Av[i,j] - end + copytrito_upper!(B, A) else # uplo == 'L' LAPACK.lacpy_size_check(size(B), (m, m < n ? m : n)) - for j in axes(A,2), i in axes(A,1)[j:end] - # extract the parents for LowerTriangular matrices - Bv, Av = lowertridata(B), lowertridata(A) - @inbounds Bv[i,j] = Av[i,j] - end + copytrito_lower!(B, A) end return B end +@inline function copytrito_upper!(B, A) + # extract the parents for UpperTriangular matrices + Bv, Av = uppertridata(B), uppertridata(A) + for j in axes(A,2), i in axes(A,1)[begin : min(j,end)] + @inbounds Bv[i,j] = Av[i,j] + end + B +end +@inline function copytrito_lower!(B, A) + # extract the parents for LowerTriangular matrices + Bv, Av = lowertridata(B), lowertridata(A) + for j in axes(A,2), i in axes(A,1)[j:end] + @inbounds Bv[i,j] = Av[i,j] + end + B +end # Forward LAPACK-compatible strided matrices to lacpy function copytrito!(B::StridedMatrixStride1{T}, A::StridedMatrixStride1{T}, uplo::AbstractChar) where {T<:BlasFloat} LAPACK.lacpy!(B, A, uplo) diff --git a/src/triangular.jl b/src/triangular.jl index ad66f492..9149c3d4 100644 --- a/src/triangular.jl +++ b/src/triangular.jl @@ -152,6 +152,12 @@ lowertriangular(U::LowerOrUnitLowerTriangular) = U Base.dataids(A::UpperOrLowerTriangular) = Base.dataids(A.data) +function Matrix{T}(U::UpperOrLowerTriangular) where {T} + M = Matrix{T}(undef, size(U)) + copyto_unaliased!(M, U) + return M +end + imag(A::UpperTriangular) = UpperTriangular(imag(A.data)) imag(A::LowerTriangular) = LowerTriangular(imag(A.data)) imag(A::UpperTriangular{<:Any,<:StridedMaybeAdjOrTransMat}) = imag.(A) From aa9bf03f866fd2d4a55e01cf8de8861b11cb9bec Mon Sep 17 00:00:00 2001 From: Jishnu Bhattacharya Date: Thu, 10 Apr 2025 16:34:07 +0530 Subject: [PATCH 2/8] Remove branches in copyto_unaliased --- src/triangular.jl | 44 ++++++++++++++++++++++++++++++++------------ 1 file changed, 32 insertions(+), 12 deletions(-) diff --git a/src/triangular.jl b/src/triangular.jl index 9149c3d4..9d3b8f13 100644 --- a/src/triangular.jl +++ b/src/triangular.jl @@ -628,29 +628,49 @@ end end # for strided matrices, we explicitly loop over the arrays to improve cache locality # This fuses the copytrito! for the two halves +@propagate_inbounds function copy_unaliased_stored!(dest, U::UpperTriangular, col) + for row in firstindex(dest,1):col + dest[row,col] = U.data[row,col] + end + dest +end +@propagate_inbounds function copy_unaliased_stored!(dest, U::UnitUpperTriangular, col) + for row in firstindex(dest,1):col-1 + dest[row,col] = U.data[row,col] + end + dest[col, col] = U[BandIndex(0,col)] + dest +end +@propagate_inbounds function copy_unaliased_stored!(dest, U::LowerTriangular, col) + for row in col:lastindex(dest,1) + dest[row,col] = L.data[row,col] + end + dest +end +@propagate_inbounds function copy_unaliased_stored!(dest, U::UnitLowerTriangular, col) + dest[col, col] = U[BandIndex(0,col)] + for row in col+1:lastindex(dest,1) + dest[row,col] = L.data[row,col] + end + dest +end @inline function copy_unaliased!(dest::StridedMatrix, U::UpperOrUnitUpperTriangular{<:Any, <:StridedMatrix}) @boundscheck checkbounds(dest, axes(U)...) - isunit = U isa UnitUpperTriangular for col in axes(dest,2) - for row in firstindex(dest,1):col-isunit - @inbounds dest[row,col] = U.data[row,col] - end - for row in col+!isunit:lastindex(dest,1) - @inbounds dest[row,col] = U[row,col] + @inbounds copy_unaliased_stored!(dest, U, col) + for row in col+1:lastindex(dest,1) + @inbounds dest[row,col] = diagzero(U,row,col) end end return dest end @inline function copy_unaliased!(dest::StridedMatrix, L::LowerOrUnitLowerTriangular{<:Any, <:StridedMatrix}) @boundscheck checkbounds(dest, axes(L)...) - isunit = L isa UnitLowerTriangular for col in axes(dest,2) - for row in firstindex(dest,1):col-!isunit - @inbounds dest[row,col] = L[row,col] - end - for row in col+isunit:lastindex(dest,1) - @inbounds dest[row,col] = L.data[row,col] + for row in firstindex(dest,1):col-1 + @inbounds dest[row,col] = diagzero(L,row,col) end + @inbounds copy_unaliased_stored!(dest, L, col) end return dest end From 80062b50dc2804114b17af89c103cd4955bc0bf2 Mon Sep 17 00:00:00 2001 From: Jishnu Bhattacharya Date: Thu, 10 Apr 2025 16:46:11 +0530 Subject: [PATCH 3/8] Undo changes to `copytrito!` --- src/generic.jl | 28 ++++++++++------------------ 1 file changed, 10 insertions(+), 18 deletions(-) diff --git a/src/generic.jl b/src/generic.jl index 4dea4643..b99dc6df 100644 --- a/src/generic.jl +++ b/src/generic.jl @@ -2079,29 +2079,21 @@ function copytrito!(B::AbstractMatrix, A::AbstractMatrix, uplo::AbstractChar) A = Base.unalias(B, A) if uplo == 'U' LAPACK.lacpy_size_check(size(B), (n < m ? n : m, n)) - copytrito_upper!(B, A) + # extract the parents for UpperTriangular matrices + Bv, Av = uppertridata(B), uppertridata(A) + for j in axes(A,2), i in axes(A,1)[begin : min(j,end)] + @inbounds Bv[i,j] = Av[i,j] + end else # uplo == 'L' LAPACK.lacpy_size_check(size(B), (m, m < n ? m : n)) - copytrito_lower!(B, A) + # extract the parents for LowerTriangular matrices + Bv, Av = lowertridata(B), lowertridata(A) + for j in axes(A,2), i in axes(A,1)[j:end] + @inbounds Bv[i,j] = Av[i,j] + end end return B end -@inline function copytrito_upper!(B, A) - # extract the parents for UpperTriangular matrices - Bv, Av = uppertridata(B), uppertridata(A) - for j in axes(A,2), i in axes(A,1)[begin : min(j,end)] - @inbounds Bv[i,j] = Av[i,j] - end - B -end -@inline function copytrito_lower!(B, A) - # extract the parents for LowerTriangular matrices - Bv, Av = lowertridata(B), lowertridata(A) - for j in axes(A,2), i in axes(A,1)[j:end] - @inbounds Bv[i,j] = Av[i,j] - end - B -end # Forward LAPACK-compatible strided matrices to lacpy function copytrito!(B::StridedMatrixStride1{T}, A::StridedMatrixStride1{T}, uplo::AbstractChar) where {T<:BlasFloat} LAPACK.lacpy!(B, A, uplo) From 950959b4c68d2ba70d7628c43b5d0e86fbfc395f Mon Sep 17 00:00:00 2001 From: Jishnu Bhattacharya Date: Thu, 10 Apr 2025 17:16:53 +0530 Subject: [PATCH 4/8] fix variable names --- src/triangular.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/triangular.jl b/src/triangular.jl index 9d3b8f13..3c3b57bb 100644 --- a/src/triangular.jl +++ b/src/triangular.jl @@ -641,14 +641,14 @@ end dest[col, col] = U[BandIndex(0,col)] dest end -@propagate_inbounds function copy_unaliased_stored!(dest, U::LowerTriangular, col) +@propagate_inbounds function copy_unaliased_stored!(dest, L::LowerTriangular, col) for row in col:lastindex(dest,1) dest[row,col] = L.data[row,col] end dest end -@propagate_inbounds function copy_unaliased_stored!(dest, U::UnitLowerTriangular, col) - dest[col, col] = U[BandIndex(0,col)] +@propagate_inbounds function copy_unaliased_stored!(dest, L::UnitLowerTriangular, col) + dest[col, col] = L[BandIndex(0,col)] for row in col+1:lastindex(dest,1) dest[row,col] = L.data[row,col] end From a28a14ed48a34791509915b1dc87df3831807103 Mon Sep 17 00:00:00 2001 From: Jishnu Bhattacharya Date: Thu, 10 Apr 2025 18:03:32 +0530 Subject: [PATCH 5/8] Call `_copyto!` instead of `copyto_unaliased!` --- src/triangular.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/triangular.jl b/src/triangular.jl index 3c3b57bb..d6e05341 100644 --- a/src/triangular.jl +++ b/src/triangular.jl @@ -154,7 +154,7 @@ Base.dataids(A::UpperOrLowerTriangular) = Base.dataids(A.data) function Matrix{T}(U::UpperOrLowerTriangular) where {T} M = Matrix{T}(undef, size(U)) - copyto_unaliased!(M, U) + _copyto!(M, U) return M end From a41b019a3d45ba8e851bad195b880780d91eb242 Mon Sep 17 00:00:00 2001 From: Jishnu Bhattacharya Date: Sun, 13 Apr 2025 16:55:36 +0530 Subject: [PATCH 6/8] Use `copy!` in `Matrix` constructor --- src/triangular.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/triangular.jl b/src/triangular.jl index d6e05341..b957fd10 100644 --- a/src/triangular.jl +++ b/src/triangular.jl @@ -154,7 +154,7 @@ Base.dataids(A::UpperOrLowerTriangular) = Base.dataids(A.data) function Matrix{T}(U::UpperOrLowerTriangular) where {T} M = Matrix{T}(undef, size(U)) - _copyto!(M, U) + copy!(M, U) return M end From 72e5cc84d8dcc2849a60876c183fd6e10a3838b4 Mon Sep 17 00:00:00 2001 From: Jishnu Bhattacharya Date: Sun, 13 Apr 2025 16:57:13 +0530 Subject: [PATCH 7/8] Revert unnecessary changes --- src/generic.jl | 28 ++++++++++++++++++---------- 1 file changed, 18 insertions(+), 10 deletions(-) diff --git a/src/generic.jl b/src/generic.jl index b99dc6df..4dea4643 100644 --- a/src/generic.jl +++ b/src/generic.jl @@ -2079,21 +2079,29 @@ function copytrito!(B::AbstractMatrix, A::AbstractMatrix, uplo::AbstractChar) A = Base.unalias(B, A) if uplo == 'U' LAPACK.lacpy_size_check(size(B), (n < m ? n : m, n)) - # extract the parents for UpperTriangular matrices - Bv, Av = uppertridata(B), uppertridata(A) - for j in axes(A,2), i in axes(A,1)[begin : min(j,end)] - @inbounds Bv[i,j] = Av[i,j] - end + copytrito_upper!(B, A) else # uplo == 'L' LAPACK.lacpy_size_check(size(B), (m, m < n ? m : n)) - # extract the parents for LowerTriangular matrices - Bv, Av = lowertridata(B), lowertridata(A) - for j in axes(A,2), i in axes(A,1)[j:end] - @inbounds Bv[i,j] = Av[i,j] - end + copytrito_lower!(B, A) end return B end +@inline function copytrito_upper!(B, A) + # extract the parents for UpperTriangular matrices + Bv, Av = uppertridata(B), uppertridata(A) + for j in axes(A,2), i in axes(A,1)[begin : min(j,end)] + @inbounds Bv[i,j] = Av[i,j] + end + B +end +@inline function copytrito_lower!(B, A) + # extract the parents for LowerTriangular matrices + Bv, Av = lowertridata(B), lowertridata(A) + for j in axes(A,2), i in axes(A,1)[j:end] + @inbounds Bv[i,j] = Av[i,j] + end + B +end # Forward LAPACK-compatible strided matrices to lacpy function copytrito!(B::StridedMatrixStride1{T}, A::StridedMatrixStride1{T}, uplo::AbstractChar) where {T<:BlasFloat} LAPACK.lacpy!(B, A, uplo) From bef96b4e457c12f80d6ca4e3655ec4dfaa80f8ce Mon Sep 17 00:00:00 2001 From: Jishnu Bhattacharya Date: Sun, 13 Apr 2025 16:58:24 +0530 Subject: [PATCH 8/8] Redo "undo changes to `copytrito!`" --- src/generic.jl | 28 ++++++++++------------------ 1 file changed, 10 insertions(+), 18 deletions(-) diff --git a/src/generic.jl b/src/generic.jl index 4dea4643..b99dc6df 100644 --- a/src/generic.jl +++ b/src/generic.jl @@ -2079,29 +2079,21 @@ function copytrito!(B::AbstractMatrix, A::AbstractMatrix, uplo::AbstractChar) A = Base.unalias(B, A) if uplo == 'U' LAPACK.lacpy_size_check(size(B), (n < m ? n : m, n)) - copytrito_upper!(B, A) + # extract the parents for UpperTriangular matrices + Bv, Av = uppertridata(B), uppertridata(A) + for j in axes(A,2), i in axes(A,1)[begin : min(j,end)] + @inbounds Bv[i,j] = Av[i,j] + end else # uplo == 'L' LAPACK.lacpy_size_check(size(B), (m, m < n ? m : n)) - copytrito_lower!(B, A) + # extract the parents for LowerTriangular matrices + Bv, Av = lowertridata(B), lowertridata(A) + for j in axes(A,2), i in axes(A,1)[j:end] + @inbounds Bv[i,j] = Av[i,j] + end end return B end -@inline function copytrito_upper!(B, A) - # extract the parents for UpperTriangular matrices - Bv, Av = uppertridata(B), uppertridata(A) - for j in axes(A,2), i in axes(A,1)[begin : min(j,end)] - @inbounds Bv[i,j] = Av[i,j] - end - B -end -@inline function copytrito_lower!(B, A) - # extract the parents for LowerTriangular matrices - Bv, Av = lowertridata(B), lowertridata(A) - for j in axes(A,2), i in axes(A,1)[j:end] - @inbounds Bv[i,j] = Av[i,j] - end - B -end # Forward LAPACK-compatible strided matrices to lacpy function copytrito!(B::StridedMatrixStride1{T}, A::StridedMatrixStride1{T}, uplo::AbstractChar) where {T<:BlasFloat} LAPACK.lacpy!(B, A, uplo)