Documentation update for LinearAlgebra functions (#52934)

Documentation update for LinearAlgebra functions, see #52725 - [X] Pivoting strategies - [X] `NoPivot` - [X] `RowNonZero` - [X] `ColumnNorm` - [X] `RowMaximum` - [X] Matrix multiplication functions - [X] `copy_transpose!` - [X] `copyto!` Note: `copyto!` is not mentioned in the original issue. However, it is felt that `copy_transpose!` without the complementing `copyto!` will be confusing for a reader. - [X] Exceptions - [X] `LAPACKException` - [X] `RankDeficientException` Co-authored-by: Steven G. Johnson <stevenj@mit.edu> Co-authored-by: Daniel Karrasch <daniel.karrasch@posteo.de>
JuliaLang · Feb 12, 2024 · f8d5947 · f8d5947
1 parent 29a58d5
commit f8d5947
Show file tree

Hide file tree

Showing 7 changed files with 188 additions and 17 deletions.
diff --git a/stdlib/LinearAlgebra/docs/src/index.md b/stdlib/LinearAlgebra/docs/src/index.md
@@ -404,6 +404,35 @@ 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 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.
+
+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.
+
+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](@ref man-linalg-factorizations) for details on the
+supported pivoting strategies.
+
+See also [`LinearAlgebra.ZeroPivotException`](@ref).
+
+```@docs
+LinearAlgebra.NoPivot
+LinearAlgebra.RowNonZero
+LinearAlgebra.RowMaximum
+LinearAlgebra.ColumnNorm
+```
+
 ## Standard functions
 
 Linear algebra functions in Julia are largely implemented by calling functions from [LAPACK](https://www.netlib.org/lapack/).
@@ -418,6 +447,8 @@ Base.:/(::AbstractVecOrMat, ::AbstractVecOrMat)
 LinearAlgebra.SingularException
 LinearAlgebra.PosDefException
 LinearAlgebra.ZeroPivotException
+LinearAlgebra.RankDeficientException
+LinearAlgebra.LAPACKException
 LinearAlgebra.dot
 LinearAlgebra.dot(::Any, ::Any, ::Any)
 LinearAlgebra.cross
@@ -571,6 +602,8 @@ LinearAlgebra.checksquare
 LinearAlgebra.peakflops
 LinearAlgebra.hermitianpart
 LinearAlgebra.hermitianpart!
+LinearAlgebra.copy_adjoint!
+LinearAlgebra.copy_transpose!
 ```
 
 ## Low-level matrix operations
@@ -762,7 +795,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/LinearAlgebra.jl b/stdlib/LinearAlgebra/src/LinearAlgebra.jl
@@ -78,6 +78,7 @@ export
     cholesky,
     cond,
     condskeel,
+    copy_adjoint!,
     copy_transpose!,
     copyto!,
     copytrito!,
@@ -190,10 +191,54 @@ abstract type Algorithm end
 struct DivideAndConquer <: Algorithm end
 struct QRIteration <: Algorithm end
 
+# 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, 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
+
+"""
+    RowNonZero
+
+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.
+"""
 struct RowNonZero <: PivotingStrategy end
+
+"""
+    RowMaximum
+
+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.
+"""
 struct RowMaximum <: PivotingStrategy end
+
+"""
+    ColumnNorm
+
+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.
+"""
 struct ColumnNorm <: PivotingStrategy end
 
 # Check that stride of matrix/vector is 1

diff --git a/stdlib/LinearAlgebra/src/exceptions.jl b/stdlib/LinearAlgebra/src/exceptions.jl
@@ -6,6 +6,13 @@ export LAPACKException,
        RankDeficientException,
        ZeroPivotException
 
+"""
+    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
@@ -41,6 +48,13 @@ function Base.showerror(io::IO, ex::PosDefException)
     print(io, "; Factorization failed.")
 end
 
+"""
+    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. The `info` field indicates the computed rank of the matrix.
+"""
 struct RankDeficientException <: Exception
     info::BlasInt
 end

diff --git a/stdlib/LinearAlgebra/src/matmul.jl b/stdlib/LinearAlgebra/src/matmul.jl
@@ -682,19 +682,60 @@ 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]` | `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) 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
 
+"""
+    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)`.
+
+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
@@ -178,8 +178,41 @@ 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)
 
-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}) -> B
+
+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)`.
+"""
+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),")")))
@@ -194,21 +227,18 @@ 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)
     end
     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(A)
+    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
@@ -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

diff --git a/stdlib/LinearAlgebra/test/runtests.jl b/stdlib/LinearAlgebra/test/runtests.jl
@@ -6,7 +6,5 @@ for file in readlines(joinpath(@__DIR__, "testgroups"))
 end
 
 @testset "Docstrings" begin
-    undoc = Docs.undocumented_names(LinearAlgebra)
-    @test_broken isempty(undoc)
-    @test undoc == [:ColumnNorm, :LAPACKException, :NoPivot, :RankDeficientException, :RowMaximum, :RowNonZero, :copy_transpose!]
+    @test isempty(Docs.undocumented_names(LinearAlgebra))
 end