Merge branch 'master' into jishnub/lrmul_adjtrans

Author: dkarrasch

docs/src/index.md

@@ -745,7 +745,7 @@ LinearAlgebra.BLAS.trsv
 
 ```@docs
 LinearAlgebra.BLAS.ger!
-# xGERU
+LinearAlgebra.BLAS.geru!
 # xGERC
 LinearAlgebra.BLAS.her!
 # xHPR

src/LinearAlgebra.jl

@@ -180,7 +180,8 @@ public AbstractTriangular,
         symmetric_type,
         zeroslike,
         matprod_dest,
-        fillstored!
+        fillstored!,
+        fillband!
 
 const BlasFloat = Union{Float64,Float32,ComplexF64,ComplexF32}
 const BlasReal = Union{Float64,Float32}

src/adjtrans.jl

@@ -580,3 +580,8 @@ diagview(A::Adjoint, k::Integer = 0) = _vecadjoint(diagview(parent(A), -k))
 # triu and tril
 triu!(A::AdjOrTransAbsMat, k::Integer = 0) = wrapperop(A)(tril!(parent(A), -k))
 tril!(A::AdjOrTransAbsMat, k::Integer = 0) = wrapperop(A)(triu!(parent(A), -k))
+
+function fillband!(A::AdjOrTrans, v, k1, k2)
+    fillband!(parent(A), wrapperop(A)(v), -k2, -k1)
+    return A
+end

src/bidiag.jl

@@ -1467,6 +1467,35 @@ function inv(B::Bidiagonal{T}) where T
     return B.uplo == 'U' ? UpperTriangular(dest) : LowerTriangular(dest)
 end
 
+# cholesky-version for (sym)tridiagonal matrices
+for (T, uplo) in ((:UpperTriangular, :(:U)), (:LowerTriangular, :(:L)))
+    @eval function _chol!(A::Bidiagonal, ::Type{$T})
+        dv = real(A.dv)
+        ev = A.ev
+        n = length(dv)
+        @inbounds for i in 1:n-1
+            iszero(dv[i]) && throw(ZeroPivotException(i))
+            ev[i] /= dv[i]
+            dv[i+1] -= abs2(ev[i])*dv[i]
+            Akk = dv[i]
+            Akk, info = _chol!(Akk, UpperTriangular)
+            if info != 0
+                return $T(A), convert(BlasInt, i)
+            end
+            dv[i] = Akk
+            ev[i] *= dv[i]
+        end
+        Akk = dv[n]
+        Akk, info = _chol!(Akk, $T)
+        if info != 0
+            return $T(A), convert(BlasInt, n)
+        end
+        dv[n] = Akk
+        B = Bidiagonal(dv, ev, $uplo)
+        return $T(B), convert(BlasInt, 0)
+    end
+end
+
 # Eigensystems
 eigvals(M::Bidiagonal) = copy(M.dv)
 function eigvecs(M::Bidiagonal{T}) where T
@@ -1543,3 +1572,30 @@ function Base._sum(A::Bidiagonal, dims::Integer)
     end
     res
 end
+
+function fillband!(B::Bidiagonal, x, l, u)
+    if l > u
+        return B
+    end
+    if ((B.uplo == 'U' && (l < 0 || u > 1)) ||
+            (B.uplo == 'L' && (l < -1 || u > 0))) && !iszero(x)
+        throw_fillband_error(l, u, x)
+    else
+        if B.uplo == 'U'
+            if l <= 1 <= u
+                fill!(B.ev, x)
+            end
+            if l <= 0 <= u
+                fill!(B.dv, x)
+            end
+        else # B.uplo == 'L'
+            if l <= 0 <= u
+                fill!(B.dv, x)
+            end
+            if l <= -1 <= u
+                fill!(B.ev, x)
+            end
+        end
+    end
+    return B
+end

src/cholesky.jl

@@ -413,7 +413,7 @@ end
 # cholesky!. Destructive methods for computing Cholesky factorization of real symmetric
 # or Hermitian matrix
 ## No pivoting (default)
-function cholesky!(A::SelfAdjoint, ::NoPivot = NoPivot(); check::Bool = true)
+function cholesky!(A::RealSymHermitian, ::NoPivot = NoPivot(); check::Bool = true)
     C, info = _chol!(A.data, A.uplo == 'U' ? UpperTriangular : LowerTriangular)
     check && checkpositivedefinite(info)
     return Cholesky(C.data, A.uplo, info)
@@ -455,7 +455,7 @@ end
 
 ## With pivoting
 ### Non BLAS/LAPACK element types (generic).
-function cholesky!(A::SelfAdjoint, ::RowMaximum; tol = 0.0, check::Bool = true)
+function cholesky!(A::RealSymHermitian, ::RowMaximum; tol = 0.0, check::Bool = true)
     AA, piv, rank, info = _cholpivoted!(A.data, A.uplo == 'U' ? UpperTriangular : LowerTriangular, tol, check)
     C = CholeskyPivoted(AA, A.uplo, piv, rank, tol, info)
     check && chkfullrank(C)

src/dense.jl

@@ -205,6 +205,23 @@ tril(M::Matrix, k::Integer) = tril!(copy(M), k)
     fillband!(A::AbstractMatrix, x, l, u)
 
 Fill the band between diagonals `l` and `u` with the value `x`.
+
+# Examples
+```jldoctest
+julia> A = zeros(4,4)
+4×4 Matrix{Float64}:
+ 0.0  0.0  0.0  0.0
+ 0.0  0.0  0.0  0.0
+ 0.0  0.0  0.0  0.0
+ 0.0  0.0  0.0  0.0
+
+julia> LinearAlgebra.fillband!(A, 2, 0, 1)
+4×4 Matrix{Float64}:
+ 2.0  2.0  0.0  0.0
+ 0.0  2.0  2.0  0.0
+ 0.0  0.0  2.0  2.0
+ 0.0  0.0  0.0  2.0
+```
 """
 function fillband!(A::AbstractMatrix{T}, x, l, u) where T
     require_one_based_indexing(A)
@@ -1731,10 +1748,21 @@ both with the value of `M` and the intended application of the pseudoinverse.
 The default relative tolerance is `n*ϵ`, where `n` is the size of the smallest
 dimension of `M`, and `ϵ` is the [`eps`](@ref) of the element type of `M`.
 
-For inverting dense ill-conditioned matrices in a least-squares sense,
-`rtol = sqrt(eps(real(float(oneunit(eltype(M))))))` is recommended.
+For solving dense, ill-conditioned equations in a least-square sense, it
+is better to *not* explicitly form the pseudoinverse matrix, since this
+can lead to numerical instability at low tolerances.  The default `M \\ b`
+algorithm instead uses pivoted QR factorization ([`qr`](@ref)).  To use an
+SVD-based algorithm, it is better to employ the SVD directly via `svd(M; rtol, atol) \\ b`
+or `ldiv!(svd(M), b; rtol, atol)`.
+
+One can also pass `M = svd(A)` as the argument to `pinv` in order to re-use
+an existing [`SVD`](@ref) factorization.  In this case, `pinv` will return
+the SVD of the pseudo-inverse, which can be applied accurately, instead of an explicit matrix.
+
+!!! compat "Julia 1.13"
+    Passing an `SVD` object to `pinv` requires Julia 1.13 or later.
 
-For more information, see [^issue8859], [^B96], [^S84], [^KY88].
+For more information, see [^pr1387], [^B96], [^S84], [^KY88].
 
 # Examples
 ```jldoctest
@@ -1754,15 +1782,15 @@ julia> M * N
  4.44089e-16   1.0
 ```
 
-[^issue8859]: Issue 8859, "Fix least squares", [https://github.com/JuliaLang/julia/pull/8859](https://github.com/JuliaLang/julia/pull/8859)
+[^pr1387]: PR 1387, "stable pinv least-squares", [LinearAlgebra.jl#1387](https://github.com/JuliaLang/LinearAlgebra.jl/pull/1387)
 
 [^B96]: Åke Björck, "Numerical Methods for Least Squares Problems",  SIAM Press, Philadelphia, 1996, "Other Titles in Applied Mathematics", Vol. 51. [doi:10.1137/1.9781611971484](http://epubs.siam.org/doi/book/10.1137/1.9781611971484)
 
 [^S84]: G. W. Stewart, "Rank Degeneracy", SIAM Journal on Scientific and Statistical Computing, 5(2), 1984, 403-413. [doi:10.1137/0905030](http://epubs.siam.org/doi/abs/10.1137/0905030)
 
 [^KY88]: Konstantinos Konstantinides and Kung Yao, "Statistical analysis of effective singular values in matrix rank determination", IEEE Transactions on Acoustics, Speech and Signal Processing, 36(5), 1988, 757-763. [doi:10.1109/29.1585](https://doi.org/10.1109/29.1585)
 """
-function pinv(A::AbstractMatrix{T}; atol::Real = 0.0, rtol::Real = (eps(real(float(oneunit(T))))*min(size(A)...))*iszero(atol)) where T
+function pinv(A::AbstractMatrix{T}; atol::Real=0, rtol::Real = (eps(real(float(oneunit(T))))*min(size(A)...))*iszero(atol)) where T
     m, n = size(A)
     Tout = typeof(zero(T)/sqrt(oneunit(T) + oneunit(T)))
     if m == 0 || n == 0
@@ -1830,7 +1858,7 @@ julia> nullspace(M, atol=0.95)
  1.0
 ```
 """
-function nullspace(A::AbstractVecOrMat; atol::Real = 0.0, rtol::Real = (min(size(A, 1), size(A, 2))*eps(real(float(oneunit(eltype(A))))))*iszero(atol))
+function nullspace(A::AbstractVecOrMat; atol::Real=0, rtol::Real = (min(size(A, 1), size(A, 2))*eps(real(float(oneunit(eltype(A))))))*iszero(atol))
     m, n = size(A, 1), size(A, 2)
     (m == 0 || n == 0) && return Matrix{eigtype(eltype(A))}(I, n, n)
     SVD = svd(A; full=true)

src/diagonal.jl

@@ -806,6 +806,11 @@ end
 
 kron(A::Diagonal, B::Diagonal) = Diagonal(kron(A.diag, B.diag))
 
+function kron!(C::Diagonal, A::Diagonal, B::Diagonal)
+    kron!(C.diag, A.diag, B.diag)
+    return C
+end
+
 function kron(A::Diagonal, B::SymTridiagonal)
     kdv = kron(A.diag, B.dv)
     # We don't need to drop the last element
@@ -1216,3 +1221,18 @@ end
 
 uppertriangular(D::Diagonal) = D
 lowertriangular(D::Diagonal) = D
+
+throw_fillband_error(l, u, x) = throw(ArgumentError(lazy"cannot set bands $l:$u to a nonzero value ($x)"))
+
+function fillband!(D::Diagonal, x, l, u)
+    if l > u
+        return D
+    end
+    if (l < 0 || u > 0) && !iszero(x)
+        throw_fillband_error(l, u, x)
+    end
+    if l <= 0 <= u
+        fill!(D.diag, x)
+    end
+    return D
+end

src/factorization.jl

@@ -110,7 +110,7 @@ Factorization{T}(A::AdjointFactorization) where {T} =
     adjoint(Factorization{T}(parent(A)))
 Factorization{T}(A::TransposeFactorization) where {T} =
     transpose(Factorization{T}(parent(A)))
-inv(F::Factorization{T}) where {T} = (n = size(F, 1); ldiv!(F, Matrix{T}(I, n, n)))
+inv(F::Factorization{T}) where {T} = (n = checksquare(F); ldiv!(F, Matrix{T}(I, n, n)))
 
 Base.hash(F::Factorization, h::UInt) = mapreduce(f -> hash(getfield(F, f)), hash, 1:nfields(F); init=h)
 Base.:(==)(  F::T, G::T) where {T<:Factorization} = all(f -> getfield(F, f) == getfield(G, f), 1:nfields(F))

src/generic.jl

@@ -1846,7 +1846,7 @@ julia> det(BigInt[1 0; 2 2]) # exact integer determinant
 function det(A::AbstractMatrix{T}) where {T}
     if istriu(A) || istril(A)
         S = promote_type(T, typeof((one(T)*zero(T) + zero(T))/one(T)))
-        return convert(S, det(UpperTriangular(A)))
+        return prod(Base.Fix1(convert, S), @view A[diagind(A)]; init=one(S))
     end
     return det(lu(A; check = false))
 end

src/hessenberg.jl

@@ -124,6 +124,17 @@ lmul!(x::Number, H::UpperHessenberg) = (lmul!(x, H.data); H)
 
 fillstored!(H::UpperHessenberg, x) = (fillband!(H.data, x, -1, size(H,2)-1); H)
 
+function fillband!(H::UpperHessenberg, x, l, u)
+    if l > u
+        return H
+    end
+    if l < -1 && !iszero(x)
+        throw_fillband_error(l, u, x)
+    end
+    fillband!(H.data, x, l, u)
+    return H
+end
+
 +(A::UpperHessenberg, B::UpperHessenberg) = UpperHessenberg(A.data+B.data)
 -(A::UpperHessenberg, B::UpperHessenberg) = UpperHessenberg(A.data-B.data)