Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

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

Merged
merged 10 commits into from
Dec 4, 2023
1 change: 1 addition & 0 deletions stdlib/LinearAlgebra/src/LinearAlgebra.jl
Original file line number Diff line number Diff line change
Expand Up @@ -76,6 +76,7 @@ export
cond,
condskeel,
copyto!,
copytrito!,
copy_transpose!,
cross,
adjoint,
Expand Down
45 changes: 45 additions & 0 deletions stdlib/LinearAlgebra/src/generic.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
53 changes: 53 additions & 0 deletions stdlib/LinearAlgebra/src/lapack.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
11 changes: 11 additions & 0 deletions stdlib/LinearAlgebra/test/generic.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
13 changes: 13 additions & 0 deletions stdlib/LinearAlgebra/test/lapack.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down