Merge branch 'master' into rank-qrpivoted-documentation

Author: Luis-Varona

src/adjtrans.jl

@@ -552,3 +552,20 @@ conj(A::Adjoint) = transpose(A.parent)
 function Base.replace_in_print_matrix(A::AdjOrTrans,i::Integer,j::Integer,s::AbstractString)
     Base.replace_in_print_matrix(parent(A), j, i, s)
 end
+
+# Special adjoint/transpose methods for Adjoint/Transpose that have been reshaped as a vector
+# these are used in transposing banded matrices by forwarding the operation to the bands
+"""
+    _vectranspose(A::AbstractVector)::AbstractVector
+
+Compute `vec(transpose(A))`, but avoid an allocating reshape if possible
+"""
+_vectranspose(A::AbstractVector) = vec(transpose(A))
+_vectranspose(A::Base.ReshapedArray{<:Any,1,<:TransposeAbsVec}) = transpose(parent(A))
+"""
+    _vecadjoint(A::AbstractVector)::AbstractVector
+
+Compute `vec(adjoint(A))`, but avoid an allocating reshape if possible
+"""
+_vecadjoint(A::AbstractVector) = vec(adjoint(A))
+_vecadjoint(A::Base.ReshapedArray{<:Any,1,<:AdjointAbsVec}) = adjoint(parent(A))

src/bidiag.jl

@@ -308,23 +308,13 @@ for func in (:conj, :copy, :real, :imag)
 end
 isreal(M::Bidiagonal) = isreal(M.dv) && isreal(M.ev)
 
-adjoint(B::Bidiagonal{<:Number}) = Bidiagonal(vec(adjoint(B.dv)), vec(adjoint(B.ev)), B.uplo == 'U' ? :L : :U)
-adjoint(B::Bidiagonal{<:Number, <:Base.ReshapedArray{<:Number,1,<:Adjoint}}) =
-    Bidiagonal(adjoint(parent(B.dv)), adjoint(parent(B.ev)), B.uplo == 'U' ? :L : :U)
-transpose(B::Bidiagonal{<:Number}) = Bidiagonal(B.dv, B.ev, B.uplo == 'U' ? :L : :U)
+adjoint(B::Bidiagonal) = Bidiagonal(_vecadjoint(B.dv), _vecadjoint(B.ev), B.uplo == 'U' ? :L : :U)
+transpose(B::Bidiagonal) = Bidiagonal(_vectranspose(B.dv), _vectranspose(B.ev), B.uplo == 'U' ? :L : :U)
 permutedims(B::Bidiagonal) = Bidiagonal(B.dv, B.ev, B.uplo == 'U' ? 'L' : 'U')
 function permutedims(B::Bidiagonal, perm)
     Base.checkdims_perm(axes(B), axes(B), perm)
     NTuple{2}(perm) == (2, 1) ? permutedims(B) : B
 end
-function Base.copy(aB::Adjoint{<:Any,<:Bidiagonal})
-    B = aB.parent
-    return Bidiagonal(map(x -> copy.(adjoint.(x)), (B.dv, B.ev))..., B.uplo == 'U' ? :L : :U)
-end
-function Base.copy(tB::Transpose{<:Any,<:Bidiagonal})
-    B = tB.parent
-    return Bidiagonal(map(x -> copy.(transpose.(x)), (B.dv, B.ev))..., B.uplo == 'U' ? :L : :U)
-end
 
 @noinline function throw_zeroband_error(A)
     uplo = A.uplo
@@ -1271,26 +1261,22 @@ function _dibimul_nonzeroalpha!(C::Bidiagonal, A::Diagonal, B::Bidiagonal, _add)
 end
 
 function mul(A::UpperOrUnitUpperTriangular, B::Bidiagonal)
-    TS = promote_op(matprod, eltype(A), eltype(B))
-    C = mul!(similar(A, TS, size(A)), A, B)
+    C = _mul(A, B)
     return B.uplo == 'U' ? UpperTriangular(C) : C
 end
 
 function mul(A::LowerOrUnitLowerTriangular, B::Bidiagonal)
-    TS = promote_op(matprod, eltype(A), eltype(B))
-    C = mul!(similar(A, TS, size(A)), A, B)
+    C = _mul(A, B)
     return B.uplo == 'L' ? LowerTriangular(C) : C
 end
 
 function mul(A::Bidiagonal, B::UpperOrUnitUpperTriangular)
-    TS = promote_op(matprod, eltype(A), eltype(B))
-    C = mul!(similar(B, TS, size(B)), A, B)
+    C = _mul(A, B)
     return A.uplo == 'U' ? UpperTriangular(C) : C
 end
 
 function mul(A::Bidiagonal, B::LowerOrUnitLowerTriangular)
-    TS = promote_op(matprod, eltype(A), eltype(B))
-    C = mul!(similar(B, TS, size(B)), A, B)
+    C = _mul(A, B)
     return A.uplo == 'L' ? LowerTriangular(C) : C
 end
 
@@ -1362,14 +1348,10 @@ function ldiv!(c::AbstractVecOrMat, A::Bidiagonal, b::AbstractVecOrMat)
     end
     return c
 end
-ldiv!(A::AdjOrTrans{<:Any,<:Bidiagonal}, b::AbstractVecOrMat) = @inline ldiv!(b, A, b)
-ldiv!(c::AbstractVecOrMat, A::AdjOrTrans{<:Any,<:Bidiagonal}, b::AbstractVecOrMat) =
-    (t = wrapperop(A); _rdiv!(t(c), t(b), t(A)); return c)
 
 ### Generic promotion methods and fallbacks
 \(A::Bidiagonal, B::AbstractVecOrMat) =
     ldiv!(matprod_dest(A, B, promote_op(\, eltype(A), eltype(B))), A, B)
-\(xA::AdjOrTrans{<:Any,<:Bidiagonal}, B::AbstractVecOrMat) = copy(xA) \ B
 
 ### Triangular specializations
 for tri in (:UpperTriangular, :UnitUpperTriangular)
@@ -1441,9 +1423,6 @@ function _rdiv!(C::AbstractMatrix, A::AbstractMatrix, B::Bidiagonal)
     C
 end
 rdiv!(A::AbstractMatrix, B::Bidiagonal) = @inline _rdiv!(A, A, B)
-rdiv!(A::AbstractMatrix, B::AdjOrTrans{<:Any,<:Bidiagonal}) = @inline _rdiv!(A, A, B)
-_rdiv!(C::AbstractMatrix, A::AbstractMatrix, B::AdjOrTrans{<:Any,<:Bidiagonal}) =
-    (t = wrapperop(B); ldiv!(t(C), t(B), t(A)); return C)
 
 /(A::AbstractMatrix, B::Bidiagonal) =
     _rdiv!(similar(A, promote_op(/, eltype(A), eltype(B)), size(A)), A, B)
@@ -1476,15 +1455,9 @@ function /(D::Diagonal, B::Bidiagonal)
     return B.uplo == 'U' ? UpperTriangular(A) : LowerTriangular(A)
 end
 
-/(A::AbstractMatrix, B::Transpose{<:Any,<:Bidiagonal}) = A / copy(B)
-/(A::AbstractMatrix, B::Adjoint{<:Any,<:Bidiagonal}) = A / copy(B)
 # disambiguation
 /(A::AdjointAbsVec, B::Bidiagonal) = adjoint(adjoint(B) \ parent(A))
 /(A::TransposeAbsVec, B::Bidiagonal) = transpose(transpose(B) \ parent(A))
-/(A::AdjointAbsVec, B::Transpose{<:Any,<:Bidiagonal}) = adjoint(adjoint(B) \ parent(A))
-/(A::TransposeAbsVec, B::Transpose{<:Any,<:Bidiagonal}) = transpose(transpose(B) \ parent(A))
-/(A::AdjointAbsVec, B::Adjoint{<:Any,<:Bidiagonal}) = adjoint(adjoint(B) \ parent(A))
-/(A::TransposeAbsVec, B::Adjoint{<:Any,<:Bidiagonal}) = transpose(transpose(B) \ parent(A))
 
 factorize(A::Bidiagonal) = A
 function inv(B::Bidiagonal{T}) where T

src/diagonal.jl

@@ -906,11 +906,8 @@ end
 end
 
 conj(D::Diagonal) = Diagonal(conj(D.diag))
-transpose(D::Diagonal{<:Number}) = D
-transpose(D::Diagonal) = Diagonal(transpose.(D.diag))
-adjoint(D::Diagonal{<:Number}) = Diagonal(vec(adjoint(D.diag)))
-adjoint(D::Diagonal{<:Number,<:Base.ReshapedArray{<:Number,1,<:Adjoint}}) = Diagonal(adjoint(parent(D.diag)))
-adjoint(D::Diagonal) = Diagonal(adjoint.(D.diag))
+transpose(D::Diagonal) = Diagonal(_vectranspose(D.diag))
+adjoint(D::Diagonal) = Diagonal(_vecadjoint(D.diag))
 permutedims(D::Diagonal) = D
 permutedims(D::Diagonal, perm) = (Base.checkdims_perm(axes(D), axes(D), perm); D)
 

src/lu.jl

@@ -652,7 +652,7 @@ end
 
     # Initialize variables
     info = 0
-    fill!(du2, 0)
+    fill!(du2, zero(eltype(du2)))
 
     @inbounds begin
         for i = 1:n

src/matmul.jl

@@ -51,11 +51,13 @@ end
 
 # Matrix-vector multiplication
 function (*)(A::StridedMaybeAdjOrTransMat{T}, x::StridedVector{S}) where {T<:BlasFloat,S<:Real}
+    matmul_size_check(size(A), size(x))
     TS = promote_op(matprod, T, S)
     y = isconcretetype(TS) ? convert(AbstractVector{TS}, x) : x
     mul!(similar(x, TS, size(A,1)), A, y)
 end
 function (*)(A::AbstractMatrix{T}, x::AbstractVector{S}) where {T,S}
+    matmul_size_check(size(A), size(x))
     TS = promote_op(matprod, T, S)
     mul!(similar(x, TS, axes(A,1)), A, x)
 end
@@ -113,7 +115,10 @@ julia> [1 1; 0 1] * [1 0; 1 1]
 """
 (*)(A::AbstractMatrix, B::AbstractMatrix) = mul(A, B)
 # we add an extra level of indirection to avoid ambiguities in *
-function mul(A::AbstractMatrix, B::AbstractMatrix)
+# We also define the core functionality within _mul to reuse the code elsewhere
+mul(A::AbstractMatrix, B::AbstractMatrix) = _mul(A, B)
+function _mul(A::AbstractMatrix, B::AbstractMatrix)
+    matmul_size_check(size(A), size(B))
     TS = promote_op(matprod, eltype(A), eltype(B))
     mul!(matprod_dest(A, B, TS), A, B)
 end

src/structuredbroadcast.jl

@@ -212,7 +212,7 @@ function copyto!(dest::Diagonal, bc::Broadcasted{<:StructuredMatrixStyle})
     isvalidstructbc(dest, bc) || return copyto!(dest, convert(Broadcasted{Nothing}, bc))
     axs = axes(dest)
     axes(bc) == axs || Broadcast.throwdm(axes(bc), axs)
-    for i in axs[1]
+    for i in eachindex(dest.diag)
         dest.diag[i] = @inbounds bc[BandIndex(0, i)]
     end
     return dest
@@ -222,15 +222,15 @@ function copyto!(dest::Bidiagonal, bc::Broadcasted{<:StructuredMatrixStyle})
     isvalidstructbc(dest, bc) || return copyto!(dest, convert(Broadcasted{Nothing}, bc))
     axs = axes(dest)
     axes(bc) == axs || Broadcast.throwdm(axes(bc), axs)
-    for i in axs[1]
+    for i in eachindex(dest.dv)
         dest.dv[i] = @inbounds bc[BandIndex(0, i)]
     end
     if dest.uplo == 'U'
-        for i = 1:size(dest, 1)-1
+        for i in eachindex(dest.ev)
             dest.ev[i] = @inbounds bc[BandIndex(1, i)]
         end
     else
-        for i = 1:size(dest, 1)-1
+        for i in eachindex(dest.ev)
             dest.ev[i] = @inbounds bc[BandIndex(-1, i)]
         end
     end
@@ -257,13 +257,13 @@ function copyto!(dest::Tridiagonal, bc::Broadcasted{<:StructuredMatrixStyle})
     isvalidstructbc(dest, bc) || return copyto!(dest, convert(Broadcasted{Nothing}, bc))
     axs = axes(dest)
     axes(bc) == axs || Broadcast.throwdm(axes(bc), axs)
-    for i in axs[1]
+    for i in eachindex(dest.d)
         dest.d[i] = @inbounds bc[BandIndex(0, i)]
     end
-    for i = 1:size(dest, 1)-1
+    for i in eachindex(dest.du)
         dest.du[i] = @inbounds bc[BandIndex(1, i)]
     end
-    for i = 1:size(dest, 1)-1
+    for i in eachindex(dest.dl)
         dest.dl[i] = @inbounds bc[BandIndex(-1, i)]
     end
     return dest

src/symmetric.jl

@@ -710,12 +710,14 @@ end
 mul(A::HermOrSym, B::HermOrSym) = A * copyto!(similar(parent(B)), B)
 # catch a few potential BLAS-cases
 function mul(A::HermOrSym{<:BlasFloat,<:StridedMatrix}, B::AdjOrTrans{<:BlasFloat,<:StridedMatrix})
+    matmul_size_check(size(A), size(B))
     T = promote_type(eltype(A), eltype(B))
     mul!(similar(B, T, (size(A, 1), size(B, 2))),
             convert(AbstractMatrix{T}, A),
             copy_oftype(B, T)) # make sure the AdjOrTrans wrapper is resolved
 end
 function mul(A::AdjOrTrans{<:BlasFloat,<:StridedMatrix}, B::HermOrSym{<:BlasFloat,<:StridedMatrix})
+    matmul_size_check(size(A), size(B))
     T = promote_type(eltype(A), eltype(B))
     mul!(similar(B, T, (size(A, 1), size(B, 2))),
             copy_oftype(A, T), # make sure the AdjOrTrans wrapper is resolved

src/symmetriceigen.jl

@@ -10,9 +10,9 @@ eigencopy_oftype(A::Symmetric{<:Complex}, S) = copyto!(similar(parent(A), S), A)
     default_eigen_alg(A)
 
 Return the default algorithm used to solve the eigensystem `A v = λ v` for a symmetric matrix `A`.
-Defaults to `LinearAlegbra.DivideAndConquer()`, which corresponds to the LAPACK function `LAPACK.syevd!`.
+Defaults to `LinearAlegbra.RobustRepresentations()`, which corresponds to the LAPACK function `LAPACK.syevr!`.
 """
-default_eigen_alg(@nospecialize(A)) = DivideAndConquer()
+default_eigen_alg(@nospecialize(A)) = RobustRepresentations()
 
 # Eigensolvers for symmetric and Hermitian matrices
 function eigen!(A::RealHermSymComplexHerm{<:BlasReal,<:StridedMatrix}; alg::Algorithm = default_eigen_alg(A), sortby::Union{Function,Nothing}=nothing)
@@ -37,9 +37,9 @@ matrix `F.vectors`. (The `k`th eigenvector can be obtained from the slice `F.vec
 Iterating the decomposition produces the components `F.values` and `F.vectors`.
 
 `alg` specifies which algorithm and LAPACK method to use for eigenvalue decomposition:
-- `alg = DivideAndConquer()` (default): Calls `LAPACK.syevd!`.
+- `alg = DivideAndConquer()`: Calls `LAPACK.syevd!`.
 - `alg = QRIteration()`: Calls `LAPACK.syev!`.
-- `alg = RobustRepresentations()`: Multiple relatively robust representations method, Calls `LAPACK.syevr!`.
+- `alg = RobustRepresentations()` (default): Multiple relatively robust representations method, Calls `LAPACK.syevr!`.
 
 See James W. Demmel et al, SIAM J. Sci. Comput. 30, 3, 1508 (2008) for
 a comparison of the accuracy and performance of different algorithms.
@@ -140,14 +140,18 @@ end
 Return the eigenvalues of `A`.
 
 `alg` specifies which algorithm and LAPACK method to use for eigenvalue decomposition:
-- `alg = DivideAndConquer()` (default): Calls `LAPACK.syevd!`.
+- `alg = DivideAndConquer()`: Calls `LAPACK.syevd!`.
 - `alg = QRIteration()`: Calls `LAPACK.syev!`.
-- `alg = RobustRepresentations()`: Multiple relatively robust representations method, Calls `LAPACK.syevr!`.
+- `alg = RobustRepresentations()` (default): Multiple relatively robust representations method, Calls `LAPACK.syevr!`.
 
 See James W. Demmel et al, SIAM J. Sci. Comput. 30, 3, 1508 (2008) for
 a comparison of the accuracy and performance of different methods.
 
 The default `alg` used may change in the future.
+
+!!! compat "Julia 1.12"
+    The `alg` keyword argument requires Julia 1.12 or later.
+
 """
 function eigvals(A::RealHermSymComplexHerm; alg::Algorithm = default_eigen_alg(A), sortby::Union{Function,Nothing}=nothing)
     S = eigtype(eltype(A))

src/uniformscaling.jl

@@ -217,6 +217,17 @@ function (-)(J::UniformScaling{<:Complex}, A::Hermitian)
     return B
 end
 
+function (+)(A::AdjOrTransAbsMat, J::UniformScaling)
+    checksquare(A)
+    op = wrapperop(A)
+    op(op(A) + op(J))
+end
+function (-)(J::UniformScaling, A::AdjOrTransAbsMat)
+    checksquare(A)
+    op = wrapperop(A)
+    op(op(J) - op(A))
+end
+
 function (+)(A::AbstractMatrix, J::UniformScaling)
     checksquare(A)
     B = copymutable_oftype(A, Base.promote_op(+, eltype(A), typeof(J)))
@@ -338,7 +349,8 @@ function ==(A::AbstractMatrix, J::UniformScaling)
     size(A, 1) == size(A, 2) || return false
     iszero(J.λ) && return iszero(A)
     isone(J.λ) && return isone(A)
-    return A == J.λ*one(A)
+    isdiag(A) || return false
+    return all(==(J.λ), diagview(A))
 end
 function ==(A::StridedMatrix, J::UniformScaling)
     size(A, 1) == size(A, 2) || return false

test/adjtrans.jl

@@ -6,12 +6,13 @@ isdefined(Main, :pruned_old_LA) || @eval Main include("prune_old_LA.jl")
 
 using Test, LinearAlgebra
 
-const BASE_TEST_PATH = joinpath(Sys.BINDIR, "..", "share", "julia", "test")
+const BASE_TEST_PATH = joinpath(dirname(pathof(LinearAlgebra)), "..", "test")
+const TESTHELPERS = joinpath(BASE_TEST_PATH, "testhelpers")
 
-isdefined(Main, :OffsetArrays) || @eval Main include(joinpath($(BASE_TEST_PATH), "testhelpers", "OffsetArrays.jl"))
+isdefined(Main, :OffsetArrays) || @eval Main include(joinpath($TESTHELPERS, "OffsetArrays.jl"))
 using .Main.OffsetArrays
 
-isdefined(Main, :ImmutableArrays) || @eval Main include(joinpath($(BASE_TEST_PATH), "testhelpers", "ImmutableArrays.jl"))
+isdefined(Main, :ImmutableArrays) || @eval Main include(joinpath($TESTHELPERS, "ImmutableArrays.jl"))
 using .Main.ImmutableArrays
 
 @testset "Adjoint and Transpose inner constructor basics" begin