From fb9bdd91a34b578e31c8c2e5812602dbdbe3e687 Mon Sep 17 00:00:00 2001 From: Aravindh Krishnamoorthy Date: Wed, 17 Jan 2024 00:02:32 +0100 Subject: [PATCH 01/18] Addnl LA Documentation: Dummy commit with TODOs to get a PR started. --- stdlib/LinearAlgebra/src/LinearAlgebra.jl | 2 ++ stdlib/LinearAlgebra/src/exceptions.jl | 4 ++++ stdlib/LinearAlgebra/src/matmul.jl | 2 ++ 3 files changed, 8 insertions(+) diff --git a/stdlib/LinearAlgebra/src/LinearAlgebra.jl b/stdlib/LinearAlgebra/src/LinearAlgebra.jl index bace428c1387c..30de3961fe2b7 100644 --- a/stdlib/LinearAlgebra/src/LinearAlgebra.jl +++ b/stdlib/LinearAlgebra/src/LinearAlgebra.jl @@ -190,6 +190,8 @@ abstract type Algorithm end struct DivideAndConquer <: Algorithm end struct QRIteration <: Algorithm end +# TODO: Insert docstrings for pivoting strategies here. + abstract type PivotingStrategy end struct NoPivot <: PivotingStrategy end struct RowNonZero <: PivotingStrategy end diff --git a/stdlib/LinearAlgebra/src/exceptions.jl b/stdlib/LinearAlgebra/src/exceptions.jl index 574decf79fc07..a20148c0c9d03 100644 --- a/stdlib/LinearAlgebra/src/exceptions.jl +++ b/stdlib/LinearAlgebra/src/exceptions.jl @@ -6,6 +6,8 @@ export LAPACKException, RankDeficientException, ZeroPivotException +# TODO: Insert docstring for LAPACKException here. + struct LAPACKException <: Exception info::BlasInt end @@ -41,6 +43,8 @@ function Base.showerror(io::IO, ex::PosDefException) print(io, "; Factorization failed.") end +# TODO: Insert docstring for RankDeficientException here. + struct RankDeficientException <: Exception info::BlasInt end diff --git a/stdlib/LinearAlgebra/src/matmul.jl b/stdlib/LinearAlgebra/src/matmul.jl index 9f163f8d394b4..5aeb1c3ab4147 100644 --- a/stdlib/LinearAlgebra/src/matmul.jl +++ b/stdlib/LinearAlgebra/src/matmul.jl @@ -692,6 +692,8 @@ function copyto!(B::AbstractVecOrMat, ir_dest::AbstractUnitRange{Int}, jr_dest:: B end +# TODO: Insert docstring for copy_transpose! here. + function copy_transpose!(B::AbstractMatrix, ir_dest::AbstractUnitRange{Int}, jr_dest::AbstractUnitRange{Int}, tM::AbstractChar, M::AbstractVecOrMat, ir_src::AbstractUnitRange{Int}, jr_src::AbstractUnitRange{Int}) if tM == 'N' LinearAlgebra.copy_transpose!(B, ir_dest, jr_dest, M, ir_src, jr_src) From aace7ac17125b991d12fdf636b27c6bb05d97365 Mon Sep 17 00:00:00 2001 From: Aravindh Krishnamoorthy Date: Fri, 19 Jan 2024 14:07:09 +0100 Subject: [PATCH 02/18] LinearAlgebra: Remove exceptions for exported but undocumented functions. --- stdlib/LinearAlgebra/test/runtests.jl | 2 -- 1 file changed, 2 deletions(-) diff --git a/stdlib/LinearAlgebra/test/runtests.jl b/stdlib/LinearAlgebra/test/runtests.jl index 55bc119756ce8..8e606a1d6deb6 100644 --- a/stdlib/LinearAlgebra/test/runtests.jl +++ b/stdlib/LinearAlgebra/test/runtests.jl @@ -7,6 +7,4 @@ end @testset "Docstrings" begin undoc = Docs.undocumented_names(LinearAlgebra) - @test_broken isempty(undoc) - @test undoc == [:ColumnNorm, :LAPACKException, :NoPivot, :RankDeficientException, :RowMaximum, :RowNonZero, :copy_transpose!] end From 8c8edbb041447685f9b8ac0af408884a1a61fa66 Mon Sep 17 00:00:00 2001 From: Aravindh Krishnamoorthy Date: Sat, 20 Jan 2024 10:25:22 +0100 Subject: [PATCH 03/18] Update LinearAlgebra test to ensure no undocumented exported symbols --- stdlib/LinearAlgebra/test/runtests.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/stdlib/LinearAlgebra/test/runtests.jl b/stdlib/LinearAlgebra/test/runtests.jl index 8e606a1d6deb6..61d1efcc5930b 100644 --- a/stdlib/LinearAlgebra/test/runtests.jl +++ b/stdlib/LinearAlgebra/test/runtests.jl @@ -7,4 +7,5 @@ end @testset "Docstrings" begin undoc = Docs.undocumented_names(LinearAlgebra) + @test undoc == [] end From 7bee3aef492f39c1831df0cf4199e567c62132a7 Mon Sep 17 00:00:00 2001 From: Aravindh Krishnamoorthy Date: Sat, 20 Jan 2024 17:17:27 +0100 Subject: [PATCH 04/18] Update stdlib/LinearAlgebra/test/runtests.jl Co-authored-by: Steven G. Johnson --- stdlib/LinearAlgebra/test/runtests.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/stdlib/LinearAlgebra/test/runtests.jl b/stdlib/LinearAlgebra/test/runtests.jl index 61d1efcc5930b..a7d0e666f8c1c 100644 --- a/stdlib/LinearAlgebra/test/runtests.jl +++ b/stdlib/LinearAlgebra/test/runtests.jl @@ -6,6 +6,6 @@ for file in readlines(joinpath(@__DIR__, "testgroups")) end @testset "Docstrings" begin - undoc = Docs.undocumented_names(LinearAlgebra) + @test isempty(Docs.undocumented_names(LinearAlgebra)) @test undoc == [] end From 0f89207158bb707149a69d1fcde05a651bee1bd2 Mon Sep 17 00:00:00 2001 From: Aravindh Krishnamoorthy Date: Sat, 20 Jan 2024 17:18:07 +0100 Subject: [PATCH 05/18] Update stdlib/LinearAlgebra/test/runtests.jl --- stdlib/LinearAlgebra/test/runtests.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/stdlib/LinearAlgebra/test/runtests.jl b/stdlib/LinearAlgebra/test/runtests.jl index a7d0e666f8c1c..d64da9899ca86 100644 --- a/stdlib/LinearAlgebra/test/runtests.jl +++ b/stdlib/LinearAlgebra/test/runtests.jl @@ -7,5 +7,4 @@ end @testset "Docstrings" begin @test isempty(Docs.undocumented_names(LinearAlgebra)) - @test undoc == [] end From 4b4587c09e10aa18c703910a7fdc00ba613851ca Mon Sep 17 00:00:00 2001 From: Aravindh Krishnamoorthy Date: Wed, 24 Jan 2024 23:20:52 +0100 Subject: [PATCH 06/18] DocStrings for copy_transpose! --- stdlib/LinearAlgebra/docs/src/index.md | 1 + stdlib/LinearAlgebra/src/matmul.jl | 19 ++++++++++++++++++- stdlib/LinearAlgebra/src/transpose.jl | 11 +++++++++++ 3 files changed, 30 insertions(+), 1 deletion(-) diff --git a/stdlib/LinearAlgebra/docs/src/index.md b/stdlib/LinearAlgebra/docs/src/index.md index 138e505c79a46..92175b9e32bf7 100644 --- a/stdlib/LinearAlgebra/docs/src/index.md +++ b/stdlib/LinearAlgebra/docs/src/index.md @@ -570,6 +570,7 @@ LinearAlgebra.checksquare LinearAlgebra.peakflops LinearAlgebra.hermitianpart LinearAlgebra.hermitianpart! +LinearAlgebra.copy_transpose! ``` ## Low-level matrix operations diff --git a/stdlib/LinearAlgebra/src/matmul.jl b/stdlib/LinearAlgebra/src/matmul.jl index 5aeb1c3ab4147..88f47825a8fb1 100644 --- a/stdlib/LinearAlgebra/src/matmul.jl +++ b/stdlib/LinearAlgebra/src/matmul.jl @@ -692,8 +692,25 @@ function copyto!(B::AbstractVecOrMat, ir_dest::AbstractUnitRange{Int}, jr_dest:: B end -# TODO: Insert docstring for copy_transpose! here. +""" + copy_transpose!(B::AbstractMatrix, ir_dest::AbstractUnitRange, + jr_dest::AbstractUnitRange, tM::AbstractChar, M::AbstractVecOrMat, + ir_src::AbstractUnitRange, jr_src::AbstractUnitRange) -> B + +Efficiently copy elements of matrix `M` to `B` conditioned on the character +parameter `tM` as follows: + +| `tM` | Destination | Source | +| --- | :--- | :--- | +| `'N'` | `B[ir_dest, jr_dest]` | `transpose(M)[jr_src, ir_src]` | +| `'T'` | `B[ir_dest, jr_dest]` | `M[jr_src, ir_src]` | +| `'C'` | `B[ir_dest, jr_dest]` | `conj(M)[jr_src, ir_src]` | +The elements `B[ir_dest, jr_dest]` are overwritten. + +Furthermore, the index range parameters must satisfy `length(ir_dest) == length(jr_src)` and +`length(jr_dest) == length(ir_src)`. +""" function copy_transpose!(B::AbstractMatrix, ir_dest::AbstractUnitRange{Int}, jr_dest::AbstractUnitRange{Int}, tM::AbstractChar, M::AbstractVecOrMat, ir_src::AbstractUnitRange{Int}, jr_src::AbstractUnitRange{Int}) if tM == 'N' LinearAlgebra.copy_transpose!(B, ir_dest, jr_dest, M, ir_src, jr_src) diff --git a/stdlib/LinearAlgebra/src/transpose.jl b/stdlib/LinearAlgebra/src/transpose.jl index afc3494fc9726..3597c18b6007a 100644 --- a/stdlib/LinearAlgebra/src/transpose.jl +++ b/stdlib/LinearAlgebra/src/transpose.jl @@ -178,6 +178,17 @@ copy(::Union{Transpose,Adjoint}) Base.copy(A::TransposeAbsMat) = transpose!(similar(A.parent, reverse(axes(A.parent))), A.parent) Base.copy(A::AdjointAbsMat) = adjoint!(similar(A.parent, reverse(axes(A.parent))), A.parent) +""" + copy_transpose!(B::AbstractVecOrMat, ir_dest::AbstractRange{Int}, + jr_dest::AbstractRange{Int}, A::AbstractVecOrMat, + ir_src::AbstractRange{Int}, jr_src::AbstractRange{Int}) -> B + +Efficiently copy elements of matrix `A` to `B` with transposition: `B[ir_dest, jr_dest] = +transpose(A)[jr_src, ir_src]`. The elements `B[ir_dest, jr_dest]` are overwritten. + +Furthermore, the index range parameters must satisfy `length(ir_dest) == length(jr_src)` and +`length(jr_dest) == length(ir_src)`. +""" function copy_transpose!(B::AbstractVecOrMat, ir_dest::AbstractRange{Int}, jr_dest::AbstractRange{Int}, A::AbstractVecOrMat, ir_src::AbstractRange{Int}, jr_src::AbstractRange{Int}) if length(ir_dest) != length(jr_src) From c0acaaddd54d0c75619b82f030be2aad7e20d730 Mon Sep 17 00:00:00 2001 From: Aravindh Krishnamoorthy Date: Fri, 26 Jan 2024 21:13:15 +0100 Subject: [PATCH 07/18] DocStrings for RankDeficientException and LAPACKException. --- stdlib/LinearAlgebra/docs/src/index.md | 4 +++- stdlib/LinearAlgebra/src/exceptions.jl | 14 ++++++++++++-- 2 files changed, 15 insertions(+), 3 deletions(-) diff --git a/stdlib/LinearAlgebra/docs/src/index.md b/stdlib/LinearAlgebra/docs/src/index.md index 92175b9e32bf7..71e72bfedd726 100644 --- a/stdlib/LinearAlgebra/docs/src/index.md +++ b/stdlib/LinearAlgebra/docs/src/index.md @@ -417,6 +417,8 @@ Base.:/(::AbstractVecOrMat, ::AbstractVecOrMat) LinearAlgebra.SingularException LinearAlgebra.PosDefException LinearAlgebra.ZeroPivotException +LinearAlgebra.RankDeficientException +LinearAlgebra.LAPACKException LinearAlgebra.dot LinearAlgebra.dot(::Any, ::Any, ::Any) LinearAlgebra.cross @@ -762,7 +764,7 @@ LinearAlgebra.BLAS.trsm! LinearAlgebra.BLAS.trsm ``` -## LAPACK functions +## [LAPACK functions](@id man-linalg-lapack-functions) `LinearAlgebra.LAPACK` provides wrappers for some of the LAPACK functions for linear algebra. Those functions that overwrite one of the input arrays have names ending in `'!'`. diff --git a/stdlib/LinearAlgebra/src/exceptions.jl b/stdlib/LinearAlgebra/src/exceptions.jl index a20148c0c9d03..6391d204e675f 100644 --- a/stdlib/LinearAlgebra/src/exceptions.jl +++ b/stdlib/LinearAlgebra/src/exceptions.jl @@ -6,8 +6,13 @@ export LAPACKException, RankDeficientException, ZeroPivotException -# TODO: Insert docstring for LAPACKException here. +""" + LAPACKException +Generic LAPACK exception thrown either during direct calls to the [LAPACK functions](@ref man-linalg-lapack-functions) +or during calls to other functions that use the LAPACK functions internally but lack specialized error handling. The `info` field +contains additional information on the underlying error and depends on the LAPACK function that was invoked. +""" struct LAPACKException <: Exception info::BlasInt end @@ -43,8 +48,13 @@ function Base.showerror(io::IO, ex::PosDefException) print(io, "; Factorization failed.") end -# TODO: Insert docstring for RankDeficientException here. +""" + RankDeficientException +Exception thrown when the input matrix is [rank deficient](https://en.wikipedia.org/wiki/Rank_(linear_algebra)). Some +linear algebra functions, such as the Cholesky decomposition, are only applicable to matrices that are not rank +deficient, i.e., matrices with full rank. The `info` field indicates the computed rank of the matrix. +""" struct RankDeficientException <: Exception info::BlasInt end From f27751d1b1611e3144e8e4e7cf0094743e7a0f74 Mon Sep 17 00:00:00 2001 From: Aravindh Krishnamoorthy Date: Sun, 28 Jan 2024 00:05:39 +0100 Subject: [PATCH 08/18] Initial text on pivoting --- stdlib/LinearAlgebra/docs/src/index.md | 27 ++++++++++++++++++++++++++ 1 file changed, 27 insertions(+) diff --git a/stdlib/LinearAlgebra/docs/src/index.md b/stdlib/LinearAlgebra/docs/src/index.md index 71e72bfedd726..907c1687347d2 100644 --- a/stdlib/LinearAlgebra/docs/src/index.md +++ b/stdlib/LinearAlgebra/docs/src/index.md @@ -404,6 +404,33 @@ generally broadcasting over elements in the matrix representation fail because t be highly inefficient. For such use cases, consider computing the matrix representation up front and cache it for future reuse. +## [Pivoting Strategies](@id man-linalg-pivoting-strategies) + +Several [matrix factorizations](@ref man-linalg-factorizations) support +[pivoting](https://en.wikipedia.org/wiki/Pivot_element), where the rows and columns +of the original matrix are permuted with the intention of improving the numerical +stability of the factorization. + +Pivoting may introduce a minor computational overhead. However, several matrix +decompositions, such as the LU decomposition, may fail without pivoting. + +Unfortunately, there's no silver bullet pivoting strategy and the right approach +depends on the problem. Nevertheless the default pivoting strategy of a matrix +factorization technique is often a good first choice. + +In the following, the pivoting strategies implemented in Julia are described. Note +that not all matrix factorizations may support them. Consult the documentation of the +respective matrix factorization for details on the supported strategies. + +# TODO: Enhance and polish above text. Include DocStrings for the pivoting strategies. +# TODO: Check if a table with available pivoting strategies per factorization is necessary. +# TODO: Check if the rook pivoting strategy (bunchkaufman) is worth mentioning. +# TODO: Add pivoting strategies in the @docs section below. + +```@docs + +``` + ## Standard functions Linear algebra functions in Julia are largely implemented by calling functions from [LAPACK](https://www.netlib.org/lapack/). From de24482d09656fdab474584943de6b2e6512488f Mon Sep 17 00:00:00 2001 From: Aravindh Krishnamoorthy Date: Sun, 28 Jan 2024 19:27:51 +0100 Subject: [PATCH 09/18] Added DocStrings to copyto! --- stdlib/LinearAlgebra/src/matmul.jl | 35 +++++++++++++++++++++------ stdlib/LinearAlgebra/src/transpose.jl | 12 ++++----- 2 files changed, 34 insertions(+), 13 deletions(-) diff --git a/stdlib/LinearAlgebra/src/matmul.jl b/stdlib/LinearAlgebra/src/matmul.jl index 88f47825a8fb1..7d2f2106b4fe4 100644 --- a/stdlib/LinearAlgebra/src/matmul.jl +++ b/stdlib/LinearAlgebra/src/matmul.jl @@ -682,6 +682,26 @@ end lapack_size(t::AbstractChar, M::AbstractVecOrMat) = (size(M, t=='N' ? 1 : 2), size(M, t=='N' ? 2 : 1)) +""" + copyto!(B::AbstractMatrix, ir_dest::AbstractUnitRange, jr_dest::AbstractUnitRange, + tM::AbstractChar, + M::AbstractVecOrMat, ir_src::AbstractUnitRange, jr_src::AbstractUnitRange) -> B + +Efficiently copy elements of matrix `M` to `B` conditioned on the character +parameter `tM` as follows: + +| `tM` | Destination | Source | +| --- | :--- | :--- | +| `'N'` | `B[ir_dest, jr_dest]` | `M[ir_src, jr_src]` | +| `'T'` | `B[ir_dest, jr_dest]` | `transpose(M)[ir_src, jr_src]` | +| `'C'` | `B[ir_dest, jr_dest]` | `conj(transpose(M))[ir_src, jr_src]` | + +The elements `B[ir_dest, jr_dest]` are overwritten. Furthermore, the index range +parameters must satisfy `length(ir_dest) == length(ir_src)` and +`length(jr_dest) == length(jr_src)`. + +See also [`copy_transpose!`](@ref). +""" function copyto!(B::AbstractVecOrMat, ir_dest::AbstractUnitRange{Int}, jr_dest::AbstractUnitRange{Int}, tM::AbstractChar, M::AbstractVecOrMat, ir_src::AbstractUnitRange{Int}, jr_src::AbstractUnitRange{Int}) if tM == 'N' copyto!(B, ir_dest, jr_dest, M, ir_src, jr_src) @@ -693,9 +713,9 @@ function copyto!(B::AbstractVecOrMat, ir_dest::AbstractUnitRange{Int}, jr_dest:: end """ - copy_transpose!(B::AbstractMatrix, ir_dest::AbstractUnitRange, - jr_dest::AbstractUnitRange, tM::AbstractChar, M::AbstractVecOrMat, - ir_src::AbstractUnitRange, jr_src::AbstractUnitRange) -> B + copy_transpose!(B::AbstractMatrix, ir_dest::AbstractUnitRange, jr_dest::AbstractUnitRange, + tM::AbstractChar, + M::AbstractVecOrMat, ir_src::AbstractUnitRange, jr_src::AbstractUnitRange) -> B Efficiently copy elements of matrix `M` to `B` conditioned on the character parameter `tM` as follows: @@ -703,13 +723,14 @@ parameter `tM` as follows: | `tM` | Destination | Source | | --- | :--- | :--- | | `'N'` | `B[ir_dest, jr_dest]` | `transpose(M)[jr_src, ir_src]` | -| `'T'` | `B[ir_dest, jr_dest]` | `M[jr_src, ir_src]` | +| `'T'` | `B[ir_dest, jr_dest]` | `M[jr_src, ir_src]` | | `'C'` | `B[ir_dest, jr_dest]` | `conj(M)[jr_src, ir_src]` | -The elements `B[ir_dest, jr_dest]` are overwritten. - -Furthermore, the index range parameters must satisfy `length(ir_dest) == length(jr_src)` and +The elements `B[ir_dest, jr_dest]` are overwritten. Furthermore, the index +range parameters must satisfy `length(ir_dest) == length(jr_src)` and `length(jr_dest) == length(ir_src)`. + +See also [`copyto!`](@ref). """ function copy_transpose!(B::AbstractMatrix, ir_dest::AbstractUnitRange{Int}, jr_dest::AbstractUnitRange{Int}, tM::AbstractChar, M::AbstractVecOrMat, ir_src::AbstractUnitRange{Int}, jr_src::AbstractUnitRange{Int}) if tM == 'N' diff --git a/stdlib/LinearAlgebra/src/transpose.jl b/stdlib/LinearAlgebra/src/transpose.jl index 3597c18b6007a..228c203aa1db6 100644 --- a/stdlib/LinearAlgebra/src/transpose.jl +++ b/stdlib/LinearAlgebra/src/transpose.jl @@ -179,15 +179,15 @@ Base.copy(A::TransposeAbsMat) = transpose!(similar(A.parent, reverse(axes(A.pare Base.copy(A::AdjointAbsMat) = adjoint!(similar(A.parent, reverse(axes(A.parent))), A.parent) """ - copy_transpose!(B::AbstractVecOrMat, ir_dest::AbstractRange{Int}, - jr_dest::AbstractRange{Int}, A::AbstractVecOrMat, - ir_src::AbstractRange{Int}, jr_src::AbstractRange{Int}) -> B + copy_transpose!(B::AbstractVecOrMat, ir_dest::AbstractRange{Int}, jr_dest::AbstractRange{Int}, + A::AbstractVecOrMat, ir_src::AbstractRange{Int}, jr_src::AbstractRange{Int}) -> B Efficiently copy elements of matrix `A` to `B` with transposition: `B[ir_dest, jr_dest] = -transpose(A)[jr_src, ir_src]`. The elements `B[ir_dest, jr_dest]` are overwritten. - -Furthermore, the index range parameters must satisfy `length(ir_dest) == length(jr_src)` and +transpose(A)[jr_src, ir_src]`. The elements `B[ir_dest, jr_dest]` are overwritten. Furthermore, +the index range parameters must satisfy `length(ir_dest) == length(jr_src)` and `length(jr_dest) == length(ir_src)`. + +See also [`copyto!`](@ref). """ function copy_transpose!(B::AbstractVecOrMat, ir_dest::AbstractRange{Int}, jr_dest::AbstractRange{Int}, A::AbstractVecOrMat, ir_src::AbstractRange{Int}, jr_src::AbstractRange{Int}) From 6f3338ee1bb4a1f89c4bb248f0acc7f0a836eb13 Mon Sep 17 00:00:00 2001 From: Aravindh Krishnamoorthy Date: Mon, 29 Jan 2024 21:09:49 +0100 Subject: [PATCH 10/18] Documentation for pivoting and pivoting strategies. --- stdlib/LinearAlgebra/docs/src/index.md | 28 +++++++-------- stdlib/LinearAlgebra/src/LinearAlgebra.jl | 43 +++++++++++++++++++++-- 2 files changed, 53 insertions(+), 18 deletions(-) diff --git a/stdlib/LinearAlgebra/docs/src/index.md b/stdlib/LinearAlgebra/docs/src/index.md index 907c1687347d2..ec908316f8e52 100644 --- a/stdlib/LinearAlgebra/docs/src/index.md +++ b/stdlib/LinearAlgebra/docs/src/index.md @@ -406,29 +406,25 @@ up front and cache it for future reuse. ## [Pivoting Strategies](@id man-linalg-pivoting-strategies) -Several [matrix factorizations](@ref man-linalg-factorizations) support -[pivoting](https://en.wikipedia.org/wiki/Pivot_element), where the rows and columns -of the original matrix are permuted with the intention of improving the numerical -stability of the factorization. +Several [matrix factorizations](@ref man-linalg-factorizations) support pivoting, +where a [pivot element](https://en.wikipedia.org/wiki/Pivot_element) is selected +by permuting the rows and columns of the original matrix according to a pivoting +strategy. Then, the next steps of the factorization are performed based on the +selected pivot element. -Pivoting may introduce a minor computational overhead. However, several matrix -decompositions, such as the LU decomposition, may fail without pivoting. - -Unfortunately, there's no silver bullet pivoting strategy and the right approach -depends on the problem. Nevertheless the default pivoting strategy of a matrix -factorization technique is often a good first choice. +Pivoting improves the numerical stability of the factorization. In fact, several +matrix decompositions, such as the LU decomposition, may fail without pivoting. +On the other hand, pivoting may introduce a minor computational overhead. In the following, the pivoting strategies implemented in Julia are described. Note that not all matrix factorizations may support them. Consult the documentation of the respective matrix factorization for details on the supported strategies. -# TODO: Enhance and polish above text. Include DocStrings for the pivoting strategies. -# TODO: Check if a table with available pivoting strategies per factorization is necessary. -# TODO: Check if the rook pivoting strategy (bunchkaufman) is worth mentioning. -# TODO: Add pivoting strategies in the @docs section below. - ```@docs - +LinearAlgebra.NoPivot +LinearAlgebra.RowNonZero +LinearAlgebra.RowMaximum +LinearAlgebra.ColumnNorm ``` ## Standard functions diff --git a/stdlib/LinearAlgebra/src/LinearAlgebra.jl b/stdlib/LinearAlgebra/src/LinearAlgebra.jl index 30de3961fe2b7..16808633bea57 100644 --- a/stdlib/LinearAlgebra/src/LinearAlgebra.jl +++ b/stdlib/LinearAlgebra/src/LinearAlgebra.jl @@ -190,12 +190,51 @@ abstract type Algorithm end struct DivideAndConquer <: Algorithm end struct QRIteration <: Algorithm end -# TODO: Insert docstrings for pivoting strategies here. - +# Pivoting strategies for matrix factorization algorithms. abstract type PivotingStrategy end + +""" + NoPivot + +Pivoting is not performed. Matrix factorizations such as the LU factorization +may fail without pivoting. +""" struct NoPivot <: PivotingStrategy end + +""" + RowNonZero + +First non-zero element in the remaining rows is chosen as the pivot element. +Furthermore, the rows of the matrix are permuted to bring the chosen pivot element +into the right position. + +Note that the [element type](@ref eltype) of the matrix must support the [`iszero`](@ref) +method. +""" struct RowNonZero <: PivotingStrategy end + +""" + RowMaximum + +The maximum element in the remaining rows is chosen as the pivot element. +Furthermore, the rows of the matrix are permuted to bring the chosen pivot element +into the right position. + +Note that the [element type](@ref eltype) of the matrix must support the [`abs`](@ref) +and [`<`](@ref) methods. +""" struct RowMaximum <: PivotingStrategy end + +""" + ColumnNorm + +The column with the maximum norm is used for choosing the pivot element. +Furthermore, the columns of the matrix are permuted to bring the chosen column +to the right position. + +Note that the [element type](@ref eltype) of the matrix must support the [`norm`](@ref), +[`abs`](@ref), and [`<`](@ref) methods. +""" struct ColumnNorm <: PivotingStrategy end # Check that stride of matrix/vector is 1 From 61bf648519edb0a7b686c42e7fc23556cfc61f82 Mon Sep 17 00:00:00 2001 From: Aravindh Krishnamoorthy Date: Sat, 3 Feb 2024 23:38:26 +0100 Subject: [PATCH 11/18] Incorporate changes from self review. --- stdlib/LinearAlgebra/docs/src/index.md | 26 ++++++++++++++--------- stdlib/LinearAlgebra/src/LinearAlgebra.jl | 8 +------ stdlib/LinearAlgebra/src/exceptions.jl | 2 +- stdlib/LinearAlgebra/src/transpose.jl | 9 ++++---- 4 files changed, 23 insertions(+), 22 deletions(-) diff --git a/stdlib/LinearAlgebra/docs/src/index.md b/stdlib/LinearAlgebra/docs/src/index.md index ec908316f8e52..84fb9809b2894 100644 --- a/stdlib/LinearAlgebra/docs/src/index.md +++ b/stdlib/LinearAlgebra/docs/src/index.md @@ -406,19 +406,25 @@ up front and cache it for future reuse. ## [Pivoting Strategies](@id man-linalg-pivoting-strategies) -Several [matrix factorizations](@ref man-linalg-factorizations) support pivoting, -where a [pivot element](https://en.wikipedia.org/wiki/Pivot_element) is selected -by permuting the rows and columns of the original matrix according to a pivoting -strategy. Then, the next steps of the factorization are performed based on the -selected pivot element. +Several of Julia's [matrix factorizations](@ref man-linalg-factorizations) support +[pivoting](https://en.wikipedia.org/wiki/Pivot_element), which can be used to improve their +numerical stability. In fact, some matrix factorizations, such as the LU +factorization, may fail without pivoting. -Pivoting improves the numerical stability of the factorization. In fact, several -matrix decompositions, such as the LU decomposition, may fail without pivoting. -On the other hand, pivoting may introduce a minor computational overhead. +In pivoting, first, a [pivot element](https://en.wikipedia.org/wiki/Pivot_element) +with good numerical properties is chosen based on a pivoting strategy. Next, the rows and +columns of the original matrix are permuted to bring the chosen element in place for +subsequent computation. Furthermore, the process is repeated for each stage of the factorization. -In the following, the pivoting strategies implemented in Julia are described. Note +Consequently, besides the conventional matrix factors, the outputs of +pivoted factorization schemes also include permutation matrices. + +In the following, the pivoting strategies implemented in Julia are briefly described. Note that not all matrix factorizations may support them. Consult the documentation of the -respective matrix factorization for details on the supported strategies. +respective [matrix factorization](@ref man-linalg-factorizations) for details on the +supported pivoting strategies. + +See also [`LinearAlgebra.ZeroPivotException`](@ref). ```@docs LinearAlgebra.NoPivot diff --git a/stdlib/LinearAlgebra/src/LinearAlgebra.jl b/stdlib/LinearAlgebra/src/LinearAlgebra.jl index 16808633bea57..810b330474049 100644 --- a/stdlib/LinearAlgebra/src/LinearAlgebra.jl +++ b/stdlib/LinearAlgebra/src/LinearAlgebra.jl @@ -205,8 +205,6 @@ struct NoPivot <: PivotingStrategy end RowNonZero First non-zero element in the remaining rows is chosen as the pivot element. -Furthermore, the rows of the matrix are permuted to bring the chosen pivot element -into the right position. Note that the [element type](@ref eltype) of the matrix must support the [`iszero`](@ref) method. @@ -217,8 +215,6 @@ struct RowNonZero <: PivotingStrategy end RowMaximum The maximum element in the remaining rows is chosen as the pivot element. -Furthermore, the rows of the matrix are permuted to bring the chosen pivot element -into the right position. Note that the [element type](@ref eltype) of the matrix must support the [`abs`](@ref) and [`<`](@ref) methods. @@ -228,9 +224,7 @@ struct RowMaximum <: PivotingStrategy end """ ColumnNorm -The column with the maximum norm is used for choosing the pivot element. -Furthermore, the columns of the matrix are permuted to bring the chosen column -to the right position. +The column with the maximum norm is used for subsequent computation. Note that the [element type](@ref eltype) of the matrix must support the [`norm`](@ref), [`abs`](@ref), and [`<`](@ref) methods. diff --git a/stdlib/LinearAlgebra/src/exceptions.jl b/stdlib/LinearAlgebra/src/exceptions.jl index 6391d204e675f..7791b1ddef416 100644 --- a/stdlib/LinearAlgebra/src/exceptions.jl +++ b/stdlib/LinearAlgebra/src/exceptions.jl @@ -53,7 +53,7 @@ end Exception thrown when the input matrix is [rank deficient](https://en.wikipedia.org/wiki/Rank_(linear_algebra)). Some linear algebra functions, such as the Cholesky decomposition, are only applicable to matrices that are not rank -deficient, i.e., matrices with full rank. The `info` field indicates the computed rank of the matrix. +deficient. The `info` field indicates the computed rank of the matrix. """ struct RankDeficientException <: Exception info::BlasInt diff --git a/stdlib/LinearAlgebra/src/transpose.jl b/stdlib/LinearAlgebra/src/transpose.jl index 228c203aa1db6..db2ef1574ca84 100644 --- a/stdlib/LinearAlgebra/src/transpose.jl +++ b/stdlib/LinearAlgebra/src/transpose.jl @@ -182,12 +182,13 @@ Base.copy(A::AdjointAbsMat) = adjoint!(similar(A.parent, reverse(axes(A.parent)) copy_transpose!(B::AbstractVecOrMat, ir_dest::AbstractRange{Int}, jr_dest::AbstractRange{Int}, A::AbstractVecOrMat, ir_src::AbstractRange{Int}, jr_src::AbstractRange{Int}) -> B -Efficiently copy elements of matrix `A` to `B` with transposition: `B[ir_dest, jr_dest] = -transpose(A)[jr_src, ir_src]`. The elements `B[ir_dest, jr_dest]` are overwritten. Furthermore, +Efficiently copy elements of matrix `A` to `B` with transposition as follows: + + B[ir_dest, jr_dest] = transpose(A)[jr_src, ir_src] + +The elements `B[ir_dest, jr_dest]` are overwritten. Furthermore, the index range parameters must satisfy `length(ir_dest) == length(jr_src)` and `length(jr_dest) == length(ir_src)`. - -See also [`copyto!`](@ref). """ function copy_transpose!(B::AbstractVecOrMat, ir_dest::AbstractRange{Int}, jr_dest::AbstractRange{Int}, A::AbstractVecOrMat, ir_src::AbstractRange{Int}, jr_src::AbstractRange{Int}) From d47f20e7f615c2bc3c34534f45010b20b6ac9d51 Mon Sep 17 00:00:00 2001 From: Aravindh Krishnamoorthy Date: Sat, 3 Feb 2024 23:52:17 +0100 Subject: [PATCH 12/18] index.md: Remove trailing whitespace. --- stdlib/LinearAlgebra/docs/src/index.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/stdlib/LinearAlgebra/docs/src/index.md b/stdlib/LinearAlgebra/docs/src/index.md index 84fb9809b2894..6fd584ebae28a 100644 --- a/stdlib/LinearAlgebra/docs/src/index.md +++ b/stdlib/LinearAlgebra/docs/src/index.md @@ -411,7 +411,7 @@ Several of Julia's [matrix factorizations](@ref man-linalg-factorizations) suppo numerical stability. In fact, some matrix factorizations, such as the LU factorization, may fail without pivoting. -In pivoting, first, a [pivot element](https://en.wikipedia.org/wiki/Pivot_element) +In pivoting, first, a [pivot element](https://en.wikipedia.org/wiki/Pivot_element) with good numerical properties is chosen based on a pivoting strategy. Next, the rows and columns of the original matrix are permuted to bring the chosen element in place for subsequent computation. Furthermore, the process is repeated for each stage of the factorization. From 72a188af86e54c61bbf8b5d9c59a638192be7bce Mon Sep 17 00:00:00 2001 From: Aravindh Krishnamoorthy Date: Wed, 7 Feb 2024 13:57:23 +0100 Subject: [PATCH 13/18] Apply suggestions from code review by @dkarrasch Co-authored-by: Daniel Karrasch --- stdlib/LinearAlgebra/src/LinearAlgebra.jl | 10 +++++----- stdlib/LinearAlgebra/src/matmul.jl | 2 +- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/stdlib/LinearAlgebra/src/LinearAlgebra.jl b/stdlib/LinearAlgebra/src/LinearAlgebra.jl index 810b330474049..03384e5992744 100644 --- a/stdlib/LinearAlgebra/src/LinearAlgebra.jl +++ b/stdlib/LinearAlgebra/src/LinearAlgebra.jl @@ -206,7 +206,7 @@ struct NoPivot <: PivotingStrategy end First non-zero element in the remaining rows is chosen as the pivot element. -Note that the [element type](@ref eltype) of the matrix must support the [`iszero`](@ref) +Note that the [element type](@ref eltype) of the matrix must admit an [`iszero`](@ref) method. """ struct RowNonZero <: PivotingStrategy end @@ -216,8 +216,8 @@ struct RowNonZero <: PivotingStrategy end The maximum element in the remaining rows is chosen as the pivot element. -Note that the [element type](@ref eltype) of the matrix must support the [`abs`](@ref) -and [`<`](@ref) methods. +Note that the [element type](@ref eltype) of the matrix must admit an [`abs`](@ref) method, +whose result type must admit a [`<`](@ref) method. """ struct RowMaximum <: PivotingStrategy end @@ -226,8 +226,8 @@ struct RowMaximum <: PivotingStrategy end The column with the maximum norm is used for subsequent computation. -Note that the [element type](@ref eltype) of the matrix must support the [`norm`](@ref), -[`abs`](@ref), and [`<`](@ref) methods. +Note that the [element type](@ref eltype) of the matrix must admit [`norm`](@ref) and +[`abs`](@ref) methods, whose respective result types must admit a [`<`](@ref) method. """ struct ColumnNorm <: PivotingStrategy end diff --git a/stdlib/LinearAlgebra/src/matmul.jl b/stdlib/LinearAlgebra/src/matmul.jl index 7d2f2106b4fe4..9558d2ef569db 100644 --- a/stdlib/LinearAlgebra/src/matmul.jl +++ b/stdlib/LinearAlgebra/src/matmul.jl @@ -694,7 +694,7 @@ parameter `tM` as follows: | --- | :--- | :--- | | `'N'` | `B[ir_dest, jr_dest]` | `M[ir_src, jr_src]` | | `'T'` | `B[ir_dest, jr_dest]` | `transpose(M)[ir_src, jr_src]` | -| `'C'` | `B[ir_dest, jr_dest]` | `conj(transpose(M))[ir_src, jr_src]` | +| `'C'` | `B[ir_dest, jr_dest]` | `adjoint(M)[ir_src, jr_src]` | The elements `B[ir_dest, jr_dest]` are overwritten. Furthermore, the index range parameters must satisfy `length(ir_dest) == length(ir_src)` and From e3985ddb04c29f1f4e8e82134a0e51202ec1a572 Mon Sep 17 00:00:00 2001 From: Aravindh Krishnamoorthy Date: Wed, 7 Feb 2024 20:39:50 +0100 Subject: [PATCH 14/18] Apply suggestions from code review by @stevengj and @dkarrasch Co-authored-by: Steven G. Johnson Co-authored-by: Daniel Karrasch --- stdlib/LinearAlgebra/src/LinearAlgebra.jl | 15 ++++++++++++--- stdlib/LinearAlgebra/src/matmul.jl | 2 +- 2 files changed, 13 insertions(+), 4 deletions(-) diff --git a/stdlib/LinearAlgebra/src/LinearAlgebra.jl b/stdlib/LinearAlgebra/src/LinearAlgebra.jl index 03384e5992744..3f7ac84c519db 100644 --- a/stdlib/LinearAlgebra/src/LinearAlgebra.jl +++ b/stdlib/LinearAlgebra/src/LinearAlgebra.jl @@ -197,7 +197,8 @@ abstract type PivotingStrategy end NoPivot Pivoting is not performed. Matrix factorizations such as the LU factorization -may fail without pivoting. +may fail without pivoting, and may also be numerically unstable for floating-point matrices in the face of roundoff error. +This pivot strategy is mainly useful for pedagogical purposes. """ struct NoPivot <: PivotingStrategy end @@ -206,6 +207,11 @@ struct NoPivot <: PivotingStrategy end First non-zero element in the remaining rows is chosen as the pivot element. +Beware that for floating-point matrices, the resulting LU algorithm is numerically unstable — this strategy +is mainly useful for comparison to hand calculations (which typically use this strategy) or for other +algebraic types (e.g. rational numbers) not susceptible to roundoff errors. Otherwise, the default +`RowMaximum` pivoting strategy should be generally preferred in Gaussian elimination. + Note that the [element type](@ref eltype) of the matrix must admit an [`iszero`](@ref) method. """ @@ -214,7 +220,9 @@ struct RowNonZero <: PivotingStrategy end """ RowMaximum -The maximum element in the remaining rows is chosen as the pivot element. +The maximum-magnitude element in the remaining rows is chosen as the pivot element. +This is the default strategy for LU factorization of floating-point matrices, and is sometimes +referred to as the "partial pivoting" algorithm. Note that the [element type](@ref eltype) of the matrix must admit an [`abs`](@ref) method, whose result type must admit a [`<`](@ref) method. @@ -224,7 +232,8 @@ struct RowMaximum <: PivotingStrategy end """ ColumnNorm -The column with the maximum norm is used for subsequent computation. +The column with the maximum norm is used for subsequent computation. This +is used for pivoted QR factorization. Note that the [element type](@ref eltype) of the matrix must admit [`norm`](@ref) and [`abs`](@ref) methods, whose respective result types must admit a [`<`](@ref) method. diff --git a/stdlib/LinearAlgebra/src/matmul.jl b/stdlib/LinearAlgebra/src/matmul.jl index 9558d2ef569db..7b7f25c926db2 100644 --- a/stdlib/LinearAlgebra/src/matmul.jl +++ b/stdlib/LinearAlgebra/src/matmul.jl @@ -694,7 +694,7 @@ parameter `tM` as follows: | --- | :--- | :--- | | `'N'` | `B[ir_dest, jr_dest]` | `M[ir_src, jr_src]` | | `'T'` | `B[ir_dest, jr_dest]` | `transpose(M)[ir_src, jr_src]` | -| `'C'` | `B[ir_dest, jr_dest]` | `adjoint(M)[ir_src, jr_src]` | +| `'C'` | `B[ir_dest, jr_dest]` | `conj(transpose((M))[ir_src, jr_src]` | The elements `B[ir_dest, jr_dest]` are overwritten. Furthermore, the index range parameters must satisfy `length(ir_dest) == length(ir_src)` and From 296ecd50d90e13da494d98a5512d8cf91e392513 Mon Sep 17 00:00:00 2001 From: Aravindh Krishnamoorthy Date: Wed, 7 Feb 2024 20:41:06 +0100 Subject: [PATCH 15/18] Update stdlib/LinearAlgebra/src/matmul.jl --- stdlib/LinearAlgebra/src/matmul.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/stdlib/LinearAlgebra/src/matmul.jl b/stdlib/LinearAlgebra/src/matmul.jl index 7b7f25c926db2..7d2f2106b4fe4 100644 --- a/stdlib/LinearAlgebra/src/matmul.jl +++ b/stdlib/LinearAlgebra/src/matmul.jl @@ -694,7 +694,7 @@ parameter `tM` as follows: | --- | :--- | :--- | | `'N'` | `B[ir_dest, jr_dest]` | `M[ir_src, jr_src]` | | `'T'` | `B[ir_dest, jr_dest]` | `transpose(M)[ir_src, jr_src]` | -| `'C'` | `B[ir_dest, jr_dest]` | `conj(transpose((M))[ir_src, jr_src]` | +| `'C'` | `B[ir_dest, jr_dest]` | `conj(transpose(M))[ir_src, jr_src]` | The elements `B[ir_dest, jr_dest]` are overwritten. Furthermore, the index range parameters must satisfy `length(ir_dest) == length(ir_src)` and From a6539633de2342eb8c458d6774c240b830a04be1 Mon Sep 17 00:00:00 2001 From: Aravindh Krishnamoorthy Date: Wed, 7 Feb 2024 23:29:25 +0100 Subject: [PATCH 16/18] Remove non-breaking space --- stdlib/LinearAlgebra/src/LinearAlgebra.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/stdlib/LinearAlgebra/src/LinearAlgebra.jl b/stdlib/LinearAlgebra/src/LinearAlgebra.jl index 3f7ac84c519db..7b802337fdafb 100644 --- a/stdlib/LinearAlgebra/src/LinearAlgebra.jl +++ b/stdlib/LinearAlgebra/src/LinearAlgebra.jl @@ -207,7 +207,7 @@ struct NoPivot <: PivotingStrategy end First non-zero element in the remaining rows is chosen as the pivot element. -Beware that for floating-point matrices, the resulting LU algorithm is numerically unstable — this strategy +Beware that for floating-point matrices, the resulting LU algorithm is numerically unstable — this strategy is mainly useful for comparison to hand calculations (which typically use this strategy) or for other algebraic types (e.g. rational numbers) not susceptible to roundoff errors. Otherwise, the default `RowMaximum` pivoting strategy should be generally preferred in Gaussian elimination. From 3cb092e91f120e254627b2f534aa05f8910faea2 Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Sun, 11 Feb 2024 20:25:22 +0100 Subject: [PATCH 17/18] fix adjoint case in copyto! --- stdlib/LinearAlgebra/docs/src/index.md | 1 + stdlib/LinearAlgebra/src/LinearAlgebra.jl | 1 + stdlib/LinearAlgebra/src/matmul.jl | 13 ++++---- stdlib/LinearAlgebra/src/transpose.jl | 38 +++++++++++++++++------ stdlib/LinearAlgebra/test/matmul.jl | 10 ++++++ 5 files changed, 47 insertions(+), 16 deletions(-) diff --git a/stdlib/LinearAlgebra/docs/src/index.md b/stdlib/LinearAlgebra/docs/src/index.md index 6fd584ebae28a..b8003b5b3770d 100644 --- a/stdlib/LinearAlgebra/docs/src/index.md +++ b/stdlib/LinearAlgebra/docs/src/index.md @@ -601,6 +601,7 @@ LinearAlgebra.checksquare LinearAlgebra.peakflops LinearAlgebra.hermitianpart LinearAlgebra.hermitianpart! +LinearAlgebra.copy_adjoint! LinearAlgebra.copy_transpose! ``` diff --git a/stdlib/LinearAlgebra/src/LinearAlgebra.jl b/stdlib/LinearAlgebra/src/LinearAlgebra.jl index 7b802337fdafb..8432a4d679496 100644 --- a/stdlib/LinearAlgebra/src/LinearAlgebra.jl +++ b/stdlib/LinearAlgebra/src/LinearAlgebra.jl @@ -78,6 +78,7 @@ export cholesky, cond, condskeel, + copy_adjoint!, copy_transpose!, copyto!, copytrito!, diff --git a/stdlib/LinearAlgebra/src/matmul.jl b/stdlib/LinearAlgebra/src/matmul.jl index 7d2f2106b4fe4..af676eacb046d 100644 --- a/stdlib/LinearAlgebra/src/matmul.jl +++ b/stdlib/LinearAlgebra/src/matmul.jl @@ -694,20 +694,21 @@ parameter `tM` as follows: | --- | :--- | :--- | | `'N'` | `B[ir_dest, jr_dest]` | `M[ir_src, jr_src]` | | `'T'` | `B[ir_dest, jr_dest]` | `transpose(M)[ir_src, jr_src]` | -| `'C'` | `B[ir_dest, jr_dest]` | `conj(transpose(M))[ir_src, jr_src]` | +| `'C'` | `B[ir_dest, jr_dest]` | `adjoint(M)[ir_src, jr_src]` | The elements `B[ir_dest, jr_dest]` are overwritten. Furthermore, the index range parameters must satisfy `length(ir_dest) == length(ir_src)` and `length(jr_dest) == length(jr_src)`. -See also [`copy_transpose!`](@ref). +See also [`copy_transpose!`](@ref) and [`copy_adjoint!`](@ref). """ function copyto!(B::AbstractVecOrMat, ir_dest::AbstractUnitRange{Int}, jr_dest::AbstractUnitRange{Int}, tM::AbstractChar, M::AbstractVecOrMat, ir_src::AbstractUnitRange{Int}, jr_src::AbstractUnitRange{Int}) if tM == 'N' copyto!(B, ir_dest, jr_dest, M, ir_src, jr_src) + elseif tM == 'T' + copy_transpose!(B, ir_dest, jr_dest, M, jr_src, ir_src) else - LinearAlgebra.copy_transpose!(B, ir_dest, jr_dest, M, jr_src, ir_src) - tM == 'C' && conj!(@view B[ir_dest, jr_dest]) + copy_adjoint!(B, ir_dest, jr_dest, M, jr_src, ir_src) end B end @@ -730,11 +731,11 @@ The elements `B[ir_dest, jr_dest]` are overwritten. Furthermore, the index range parameters must satisfy `length(ir_dest) == length(jr_src)` and `length(jr_dest) == length(ir_src)`. -See also [`copyto!`](@ref). +See also [`copyto!`](@ref) and [`copy_adjoint!`](@ref). """ function copy_transpose!(B::AbstractMatrix, ir_dest::AbstractUnitRange{Int}, jr_dest::AbstractUnitRange{Int}, tM::AbstractChar, M::AbstractVecOrMat, ir_src::AbstractUnitRange{Int}, jr_src::AbstractUnitRange{Int}) if tM == 'N' - LinearAlgebra.copy_transpose!(B, ir_dest, jr_dest, M, ir_src, jr_src) + copy_transpose!(B, ir_dest, jr_dest, M, ir_src, jr_src) else copyto!(B, ir_dest, jr_dest, M, jr_src, ir_src) tM == 'C' && conj!(@view B[ir_dest, jr_dest]) diff --git a/stdlib/LinearAlgebra/src/transpose.jl b/stdlib/LinearAlgebra/src/transpose.jl index db2ef1574ca84..9c1583e2cb65d 100644 --- a/stdlib/LinearAlgebra/src/transpose.jl +++ b/stdlib/LinearAlgebra/src/transpose.jl @@ -190,8 +190,29 @@ The elements `B[ir_dest, jr_dest]` are overwritten. Furthermore, the index range parameters must satisfy `length(ir_dest) == length(jr_src)` and `length(jr_dest) == length(ir_src)`. """ -function copy_transpose!(B::AbstractVecOrMat, ir_dest::AbstractRange{Int}, jr_dest::AbstractRange{Int}, - A::AbstractVecOrMat, ir_src::AbstractRange{Int}, jr_src::AbstractRange{Int}) +copy_transpose!(B::AbstractVecOrMat, ir_dest::AbstractRange{Int}, jr_dest::AbstractRange{Int}, + A::AbstractVecOrMat, ir_src::AbstractRange{Int}, jr_src::AbstractRange{Int}) = + _copy_adjtrans!(B, ir_dest, jr_dest, A, ir_src, jr_src, transpose) + +""" + copy_adjoint!(B::AbstractVecOrMat, ir_dest::AbstractRange{Int}, jr_dest::AbstractRange{Int}, + A::AbstractVecOrMat, ir_src::AbstractRange{Int}, jr_src::AbstractRange{Int}) -> B + +Efficiently copy elements of matrix `A` to `B` with adjunction as follows: + + B[ir_dest, jr_dest] = adjoint(A)[jr_src, ir_src] + +The elements `B[ir_dest, jr_dest]` are overwritten. Furthermore, +the index range parameters must satisfy `length(ir_dest) == length(jr_src)` and +`length(jr_dest) == length(ir_src)`. +""" +copy_adjoint!(B::AbstractVecOrMat, ir_dest::AbstractRange{Int}, jr_dest::AbstractRange{Int}, + A::AbstractVecOrMat, ir_src::AbstractRange{Int}, jr_src::AbstractRange{Int}) = + _copy_adjtrans!(B, ir_dest, jr_dest, A, ir_src, jr_src, adjoint) + +function _copy_adjtrans!(B::AbstractVecOrMat, ir_dest::AbstractRange{Int}, jr_dest::AbstractRange{Int}, + A::AbstractVecOrMat, ir_src::AbstractRange{Int}, jr_src::AbstractRange{Int}, + tfun::T) where {T} if length(ir_dest) != length(jr_src) throw(ArgumentError(LazyString("source and destination must have same size (got ", length(jr_src)," and ",length(ir_dest),")"))) @@ -206,7 +227,7 @@ function copy_transpose!(B::AbstractVecOrMat, ir_dest::AbstractRange{Int}, jr_de for jsrc in jr_src jdest = first(jr_dest) for isrc in ir_src - B[idest,jdest] = transpose(A[isrc,jsrc]) + B[idest,jdest] = tfun(A[isrc,jsrc]) jdest += step(jr_dest) end idest += step(ir_dest) @@ -214,13 +235,10 @@ function copy_transpose!(B::AbstractVecOrMat, ir_dest::AbstractRange{Int}, jr_de return B end -function copy_similar(A::AdjointAbsMat, ::Type{T}) where {T} - C = similar(A, T, size(A)) - adjoint!(C, parent(A)) -end -function copy_similar(A::TransposeAbsMat, ::Type{T}) where {T} - C = similar(A, T, size(A)) - transpose!(C, parent(A)) +function copy_similar(A::AdjOrTransAbsMat, ::Type{T}) where {T} + Ap = parent(A) + f! = inplace_adj_or_trans(src) + return f!(similar(Ap, T, reverse(axes(Ap))), Ap) end function Base.copyto_unaliased!(deststyle::IndexStyle, dest::AbstractMatrix, srcstyle::IndexCartesian, src::AdjOrTransAbsMat) diff --git a/stdlib/LinearAlgebra/test/matmul.jl b/stdlib/LinearAlgebra/test/matmul.jl index 4b0181f079a6f..5bcbb99a314fd 100644 --- a/stdlib/LinearAlgebra/test/matmul.jl +++ b/stdlib/LinearAlgebra/test/matmul.jl @@ -873,6 +873,16 @@ end @test Xv1' * Xv3' ≈ XcXc end +@testset "copyto! for matrices of matrices" begin + A = [randn(ComplexF64, 2,3) for _ in 1:2, _ in 1:3] + for (tfun, tM) in ((identity, 'N'), (transpose, 'T'), (adjoint, 'C')) + At = copy(tfun(A)) + B = zero.(At) + copyto!(B, axes(B, 1), axes(B, 2), tM, A, axes(A, tM == 'N' ? 1 : 2), axes(A, tM == 'N' ? 2 : 1)) + @test B == At + end +end + @testset "method ambiguity" begin # Ambiguity test is run inside a clean process. # https://github.com/JuliaLang/julia/issues/28804 From a7fd509abc1c5162def03274ae926d17cc956c3c Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Mon, 12 Feb 2024 10:34:42 +0100 Subject: [PATCH 18/18] fix typo --- stdlib/LinearAlgebra/src/transpose.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/stdlib/LinearAlgebra/src/transpose.jl b/stdlib/LinearAlgebra/src/transpose.jl index 9c1583e2cb65d..8aa04f7d34b48 100644 --- a/stdlib/LinearAlgebra/src/transpose.jl +++ b/stdlib/LinearAlgebra/src/transpose.jl @@ -237,7 +237,7 @@ end function copy_similar(A::AdjOrTransAbsMat, ::Type{T}) where {T} Ap = parent(A) - f! = inplace_adj_or_trans(src) + f! = inplace_adj_or_trans(A) return f!(similar(Ap, T, reverse(axes(Ap))), Ap) end