@@ -628,29 +628,49 @@ end
628628end
629629# for strided matrices, we explicitly loop over the arrays to improve cache locality
630630# This fuses the copytrito! for the two halves
631+ @propagate_inbounds function copy_unaliased_stored! (dest, U:: UpperTriangular , col)
632+ for row in firstindex (dest,1 ): col
633+ dest[row,col] = U. data[row,col]
634+ end
635+ dest
636+ end
637+ @propagate_inbounds function copy_unaliased_stored! (dest, U:: UnitUpperTriangular , col)
638+ for row in firstindex (dest,1 ): col- 1
639+ dest[row,col] = U. data[row,col]
640+ end
641+ dest[col, col] = U[BandIndex (0 ,col)]
642+ dest
643+ end
644+ @propagate_inbounds function copy_unaliased_stored! (dest, U:: LowerTriangular , col)
645+ for row in col: lastindex (dest,1 )
646+ dest[row,col] = L. data[row,col]
647+ end
648+ dest
649+ end
650+ @propagate_inbounds function copy_unaliased_stored! (dest, U:: UnitLowerTriangular , col)
651+ dest[col, col] = U[BandIndex (0 ,col)]
652+ for row in col+ 1 : lastindex (dest,1 )
653+ dest[row,col] = L. data[row,col]
654+ end
655+ dest
656+ end
631657@inline function copy_unaliased! (dest:: StridedMatrix , U:: UpperOrUnitUpperTriangular{<:Any, <:StridedMatrix} )
632658 @boundscheck checkbounds (dest, axes (U)... )
633- isunit = U isa UnitUpperTriangular
634659 for col in axes (dest,2 )
635- for row in firstindex (dest,1 ): col- isunit
636- @inbounds dest[row,col] = U. data[row,col]
637- end
638- for row in col+ ! isunit: lastindex (dest,1 )
639- @inbounds dest[row,col] = U[row,col]
660+ @inbounds copy_unaliased_stored! (dest, U, col)
661+ for row in col+ 1 : lastindex (dest,1 )
662+ @inbounds dest[row,col] = diagzero (U,row,col)
640663 end
641664 end
642665 return dest
643666end
644667@inline function copy_unaliased! (dest:: StridedMatrix , L:: LowerOrUnitLowerTriangular{<:Any, <:StridedMatrix} )
645668 @boundscheck checkbounds (dest, axes (L)... )
646- isunit = L isa UnitLowerTriangular
647669 for col in axes (dest,2 )
648- for row in firstindex (dest,1 ): col- ! isunit
649- @inbounds dest[row,col] = L[row,col]
650- end
651- for row in col+ isunit: lastindex (dest,1 )
652- @inbounds dest[row,col] = L. data[row,col]
670+ for row in firstindex (dest,1 ): col- 1
671+ @inbounds dest[row,col] = diagzero (L,row,col)
653672 end
673+ @inbounds copy_unaliased_stored! (dest, L, col)
654674 end
655675 return dest
656676end
0 commit comments