diff --git a/stdlib/LinearAlgebra/docs/src/index.md b/stdlib/LinearAlgebra/docs/src/index.md index 1d4c1f360feb6..da756f56e12b8 100644 --- 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 index be46a700ab44c..793775992fdea 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!, @@ -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 index 574decf79fc07..7791b1ddef416 100644 --- 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 index c84a85ebc7c81..95a24b0d798ea 100644 --- 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 index afc3494fc9726..8aa04f7d34b48 100644 --- 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,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) @@ -202,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(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 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 diff --git a/stdlib/LinearAlgebra/test/runtests.jl b/stdlib/LinearAlgebra/test/runtests.jl index 55bc119756ce8..d64da9899ca86 100644 --- 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