diff --git a/stdlib/LinearAlgebra/src/LinearAlgebra.jl b/stdlib/LinearAlgebra/src/LinearAlgebra.jl index 49f7f3226bc0a..9b250e8ee61f5 100644 --- a/stdlib/LinearAlgebra/src/LinearAlgebra.jl +++ b/stdlib/LinearAlgebra/src/LinearAlgebra.jl @@ -76,6 +76,7 @@ export cond, condskeel, copyto!, + copytrito!, copy_transpose!, cross, adjoint, diff --git a/stdlib/LinearAlgebra/src/generic.jl b/stdlib/LinearAlgebra/src/generic.jl index 9ab961461e76f..929618b95ce0a 100644 --- a/stdlib/LinearAlgebra/src/generic.jl +++ b/stdlib/LinearAlgebra/src/generic.jl @@ -1897,3 +1897,48 @@ end normalize(x) = x / norm(x) normalize(x, p::Real) = x / norm(x, p) + +""" + copytrito!(B, A, uplo) -> B + +Copies a triangular part of a matrix `A` to another matrix `B`. +`uplo` specifies the part of the matrix `A` to be copied to `B`. +Set `uplo = 'L'` for the lower triangular part or `uplo = 'U' +for the upper triangular part. + +!!! compat "Julia 1.11" + `copytrito!` requires at least Julia 1.11. + +# Examples +```jldoctest +julia> A = [1 2 ; 3 4]; + +julia> B = [0 0 ; 0 0]; + +julia> copytrito!(B, A, 'L') +2×2 Matrix{Int64}: + 1 0 + 3 4 +``` +""" +function copytrito!(B::AbstractMatrix, A::AbstractMatrix, uplo::AbstractChar) + require_one_based_indexing(A, B) + BLAS.chkuplo(uplo) + m,n = size(A) + m1,n1 = size(B) + (m1 < m || n1 < n) && throw(DimensionMismatch("B of size ($m1,$n1) should have at least the same number of rows and columns than A of size ($m,$n)")) + if uplo == 'U' + for j=1:n + for i=1:min(j,m) + @inbounds B[i,j] = A[i,j] + end + end + else # uplo == 'L' + for j=1:n + for i=j:m + @inbounds B[i,j] = A[i,j] + end + end + end + return B +end diff --git a/stdlib/LinearAlgebra/src/lapack.jl b/stdlib/LinearAlgebra/src/lapack.jl index 7bb26c2b6589a..18ebe53453e44 100644 --- a/stdlib/LinearAlgebra/src/lapack.jl +++ b/stdlib/LinearAlgebra/src/lapack.jl @@ -6992,4 +6992,57 @@ Returns `X` (overwriting `C`) and `scale`. """ trsyl!(transa::AbstractChar, transb::AbstractChar, A::AbstractMatrix, B::AbstractMatrix, C::AbstractMatrix, isgn::Int=1) +for (fn, elty) in ((:dlacpy_, :Float64), + (:slacpy_, :Float32), + (:zlacpy_, :ComplexF64), + (:clacpy_, :ComplexF32)) + @eval begin + # SUBROUTINE DLACPY( UPLO, M, N, A, LDA, B, LDB ) + # .. Scalar Arguments .. + # CHARACTER UPLO + # INTEGER LDA, LDB, M, N + # .. + # .. Array Arguments .. + # DOUBLE PRECISION A( LDA, * ), B( LDB, * ) + # .. + function lacpy!(B::AbstractMatrix{$elty}, A::AbstractMatrix{$elty}, uplo::AbstractChar) + require_one_based_indexing(A, B) + chkstride1(A, B) + m,n = size(A) + m1,n1 = size(B) + (m1 < m || n1 < n) && throw(DimensionMismatch("B of size ($m1,$n1) should have at least the same number of rows and columns than A of size ($m,$n)")) + lda = max(1, stride(A, 2)) + ldb = max(1, stride(B, 2)) + ccall((@blasfunc($fn), libblastrampoline), Cvoid, + (Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt}, Ptr{$elty}, + Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt}, Clong), + uplo, m, n, A, lda, B, ldb, 1) + B + end + end +end + +""" + lacpy!(B, A, uplo) -> B + +Copies all or part of a matrix `A` to another matrix `B`. +uplo specifies the part of the matrix `A` to be copied to `B`. +Set `uplo = 'L'` for the lower triangular part, `uplo = 'U'` +for the upper triangular part, any other character for all +the matrix `A`. + +# Examples +```jldoctest +julia> A = [1. 2. ; 3. 4.]; + +julia> B = [0. 0. ; 0. 0.]; + +julia> LAPACK.lacpy!(B, A, 'U') +2×2 Matrix{Float64}: + 1.0 2.0 + 0.0 4.0 +``` +""" +lacpy!(B::AbstractMatrix, A::AbstractMatrix, uplo::AbstractChar) + end # module diff --git a/stdlib/LinearAlgebra/test/generic.jl b/stdlib/LinearAlgebra/test/generic.jl index 303e78cfd0c58..c79173ad1011a 100644 --- a/stdlib/LinearAlgebra/test/generic.jl +++ b/stdlib/LinearAlgebra/test/generic.jl @@ -636,4 +636,15 @@ end @test condskeel(A) ≈ condskeel(A, [8,8,8]) end +@testset "copytrito!" begin + n = 10 + A = rand(n, n) + for uplo in ('L', 'U') + B = zeros(n, n) + copytrito!(B, A, uplo) + C = uplo == 'L' ? tril(A) : triu(A) + @test B ≈ C + end +end + end # module TestGeneric diff --git a/stdlib/LinearAlgebra/test/lapack.jl b/stdlib/LinearAlgebra/test/lapack.jl index 141b47d1bacef..6e12c85204a78 100644 --- a/stdlib/LinearAlgebra/test/lapack.jl +++ b/stdlib/LinearAlgebra/test/lapack.jl @@ -685,6 +685,19 @@ end end end +@testset "lacpy!" begin + @testset for elty in (Float32, Float64, ComplexF32, ComplexF64) + n = 10 + A = rand(elty, n, n) + for uplo in ('L', 'U', 'N') + B = zeros(elty, n, n) + LinearAlgebra.LAPACK.lacpy!(B, A, uplo) + C = uplo == 'L' ? tril(A) : (uplo == 'U' ? triu(A) : A) + @test B ≈ C + end + end +end + @testset "Julia vs LAPACK" begin # Test our own linear algebra functionality against LAPACK @testset for elty in (Float32, Float64, ComplexF32, ComplexF64)