[LAPACK] Interface lacpy! and add a Julia version copytrito! (#51909)

The PR adds a function `lacpy!` to easily copy a triangular part of a
square or rectangular dense matrix.

stdlib/LinearAlgebra/src/LinearAlgebra.jl

@@ -76,6 +76,7 @@ export
  cond,
  condskeel,
  copyto!,
+ copytrito!,
  copy_transpose!,
  cross,
  adjoint,

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

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

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

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)