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

Conversation

amontoison
Copy link
Contributor

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

@jishnub jishnub added the linear algebra Linear algebra label Oct 28, 2023
@jishnub
Copy link
Contributor

jishnub commented Oct 28, 2023

It'll probably be better to specialise copyto! for triangular matrices instead of adding another function. The new method may call LAPACK.lacpy! internally. I wonder how the LAPACK call performs compared to the Julia implementation?

@amontoison
Copy link
Contributor Author

amontoison commented Oct 28, 2023

Triangular wrappers only work with square matrices. It's a major issue when we want to copy the L / R factor of an LQ / QR decomposition.
But you are right, we can still specialize copyto! for triangular square matrices.

For the performance of the LAPACK version, I don't think that we have a difference for the backend OpenBLAS.
We could have a gain if we use the Intel MKL backend.

@ViralBShah
Copy link
Member

Slightly Off-topic here, but we need to have some sort of a global preference for picking a non-default BLAS.

Copy link
Member

@dkarrasch dkarrasch left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I agree with @jishnub on both points: (a) do we need the LAPACK version at all, or would it suffice to have the generic one, and (b) can't we have a more Julian name for the generic function (that may or may not call LAPACK for strided matrices of BlasFloats elements? Perhaps we could target at something closer to LinearAlgebra.copytri! (it copies a triangle into the opposite triangle of the same matrix), perhaps copytotri!?

Comment on lines 1904 to 1908
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.
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm not 100% sure, but I think the arguments/variables should be wrapped in ``, like A and uplo.

3 4
```
"""
function lacpy!(uplo::AbstractChar, A::AbstractMatrix, B::AbstractMatrix)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This cries for require_one_based_indexing, or checks on axes. The following code may then be annotated with @inbounds.

@amontoison
Copy link
Contributor Author

I like the idea copytotri!.
copytotri! is a generic Julia version that doesn't call lacpy!.
I propose to keep lacpy! because it's already interfaced.

@amontoison amontoison changed the title [LAPACK] Interface lacpy! and add a generic version [LAPACK] Interface lacpy! and add a Julia version copytotri! Oct 30, 2023
function copytotri!(uplo::AbstractChar, A::AbstractMatrix, B::AbstractMatrix)
require_one_based_indexing(A, B)
LinearAlgebra.BLAS.chkuplo(uplo)
m, n = size(A)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You definitely need to check the size of B as well, otherwise (if B is smaller than A) some indices might be out-of-bounds. Judging by the LAPACK version below, you do want to allow that B is larger than A, right?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, we should allow B to be larger than A.

stdlib/LinearAlgebra/src/generic.jl Outdated Show resolved Hide resolved
stdlib/LinearAlgebra/src/generic.jl Outdated Show resolved Hide resolved
@amontoison
Copy link
Contributor Author

@dkarrasch @jishnub The errors on CI are unrelated to this PR.
Do you have additional comments?

@ViralBShah
Copy link
Member

I'm rerunning the CI. If no other concerns, we should merge.

@amontoison amontoison force-pushed the lapack_lacpy branch 2 times, most recently from 4e1b45f to 1d0c2fc Compare November 20, 2023 22:34
@dkarrasch
Copy link
Member

dkarrasch commented Nov 27, 2023

Well, naming is hard. I had no better idea, I just thought that picking the BLAS name for a generic Julia function is not quite right, but other than that I have no strong opinions. Would copytrito! perhaps be better, like "copy a given tri angle to something else"?

@jishnub
Copy link
Contributor

jishnub commented Nov 27, 2023

This seems better, as copytri! copies the triangular part onto itself, whereas copytrito! copies it to another destination. It also reflects the fact that the structured part of the source is being copied, and not that the source is being copied to a triangular destination.

@amontoison
Copy link
Contributor Author

I renamed copytotri! into copytrito! 👍

@amontoison
Copy link
Contributor Author

@ViralBShah @dkarrasch @jishnub
Do you have any additional comments? If not, I suggest merging the pull request.

@amontoison amontoison changed the title [LAPACK] Interface lacpy! and add a Julia version copytotri! [LAPACK] Interface lacpy! and add a Julia version copytrito! Dec 4, 2023
@ViralBShah ViralBShah merged commit e280387 into JuliaLang:master Dec 4, 2023
6 of 8 checks passed
@amontoison amontoison deleted the lapack_lacpy branch December 4, 2023 13:21
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
linear algebra Linear algebra
Projects
None yet
Development

Successfully merging this pull request may close these issues.

4 participants