diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 3239dba..ed5fbd5 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -20,9 +20,7 @@ jobs: fail-fast: false matrix: version: - - '1.6' - - '1' - - '^1.9.0-0' + - '1.9' os: - ubuntu-latest - macOS-latest diff --git a/Project.toml b/Project.toml index 7ea83cf..d04beb9 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "MatrixFactorizations" uuid = "a3b82374-2e81-5b9e-98ce-41277c0e4c87" -version = "1.0" +version = "2.0" [deps] ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a" @@ -9,8 +9,8 @@ Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" [compat] -ArrayLayouts = "0.8.11, 1" -julia = "1.6" +ArrayLayouts = "1.0" +julia = "1.9" [extras] Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" diff --git a/src/MatrixFactorizations.jl b/src/MatrixFactorizations.jl index b017912..fa320be 100644 --- a/src/MatrixFactorizations.jl +++ b/src/MatrixFactorizations.jl @@ -2,13 +2,15 @@ module MatrixFactorizations using Base, LinearAlgebra, ArrayLayouts import Base: axes, axes1, getproperty, iterate, tail import LinearAlgebra: BlasInt, BlasReal, BlasFloat, BlasComplex, axpy!, - copy_oftype, checksquare, adjoint, transpose, AdjOrTrans, HermOrSym + copy_oftype, checksquare, adjoint, transpose, AdjOrTrans, HermOrSym, + det, logdet, logabsdet, isposdef import LinearAlgebra.LAPACK: chkuplo, chktrans import LinearAlgebra: cholesky, cholesky!, norm, diag, eigvals!, eigvals, eigen!, eigen, qr, axpy!, ldiv!, rdiv!, mul!, lu, lu!, ldlt, ldlt!, AbstractTriangular, inv, chkstride1, kron, lmul!, rmul!, factorize, StructuredMatrixStyle, det, logabsdet, AbstractQ, _zeros, _cut_B, _ret_size, require_one_based_indexing, checksquare, - checknonsingular, ipiv2perm, copytri!, issuccess + checknonsingular, ipiv2perm, copytri!, issuccess, RealHermSymComplexHerm, + cholcopy, checkpositivedefinite, char_uplo import Base: getindex, setindex!, *, +, -, ==, <, <=, >, >=, /, ^, \, transpose, showerror, reindex, checkbounds, @propagate_inbounds @@ -27,8 +29,8 @@ import ArrayLayouts: reflector!, reflectorApply!, materialize!, @_layoutlmul, @_ layout_getindex +export ul, ul!, ql, ql!, qrunblocked, qrunblocked!, UL, QL, reversecholesky, reversecholesky!, ReverseCholesky -export ul, ul!, ql, ql!, qrunblocked, qrunblocked!, UL, QL, choleskyinv!, choleskyinv const AdjointQtype = isdefined(LinearAlgebra, :AdjointQ) ? LinearAlgebra.AdjointQ : Adjoint const AbstractQtype = AbstractQ <: AbstractMatrix ? AbstractMatrix : AbstractQ @@ -136,7 +138,7 @@ include("ul.jl") include("qr.jl") include("ql.jl") include("rq.jl") -include("choleskyinv.jl") include("polar.jl") +include("reversecholesky.jl") end #module diff --git a/src/reversecholesky.jl b/src/reversecholesky.jl new file mode 100644 index 0000000..a04802f --- /dev/null +++ b/src/reversecholesky.jl @@ -0,0 +1,309 @@ +# This file was modified from +# a part of Julia. License is MIT: https://julialang.org/license + +########################## +# Reverse Cholesky Factorization # +########################## + +""" + ReverseCholesky <: Factorization + +Matrix factorization type of the reverse Cholesky factorization of a dense symmetric/Hermitian +positive definite matrix `A`. This is the return type of [`reversecholesky`](@ref), +the corresponding matrix factorization function. + +The triangular reverse Cholesky factor can be obtained from the factorization `F::ReverseCholesky` +via `F.L` and `F.U`, where `A ≈ F.U * F.U' ≈ F.L' * F.L`. +``` +""" +struct ReverseCholesky{T,S<:AbstractMatrix} <: Factorization{T} + factors::S + uplo::Char + info::BlasInt + + function ReverseCholesky{T,S}(factors, uplo, info) where {T,S<:AbstractMatrix} + require_one_based_indexing(factors) + new(factors, uplo, info) + end +end +ReverseCholesky(A::AbstractMatrix{T}, uplo::Symbol, info::Integer) where {T} = + ReverseCholesky{T,typeof(A)}(A, char_uplo(uplo), info) +ReverseCholesky(A::AbstractMatrix{T}, uplo::AbstractChar, info::Integer) where {T} = + ReverseCholesky{T,typeof(A)}(A, uplo, info) +ReverseCholesky(U::UpperTriangular{T}) where {T} = ReverseCholesky{T,typeof(U.data)}(U.data, 'U', 0) +ReverseCholesky(L::LowerTriangular{T}) where {T} = ReverseCholesky{T,typeof(L.data)}(L.data, 'L', 0) + +# iteration for destructuring into components +Base.iterate(C::ReverseCholesky) = (C.U, Val(:U)) +Base.iterate(C::ReverseCholesky, ::Val{:U}) = (C.L, Val(:done)) +Base.iterate(C::ReverseCholesky, ::Val{:done}) = nothing + +## Non BLAS/LAPACK element types (generic) +function _reverse_chol!(A::AbstractMatrix, ::Type{UpperTriangular}) + require_one_based_indexing(A) + n = checksquare(A) + realdiag = eltype(A) <: Complex + @inbounds begin + for k = n:-1:1 + cs = colsupport(A, k) + rs = rowsupport(A, k) + Akk = realdiag ? real(A[k,k]) : A[k,k] + for j = (k+1:n) ∩ rs + Akk -= realdiag ? abs2(A[k,j]) : A[k,j]'A[k,j] + end + A[k,k] = Akk + Akk, info = _reverse_chol!(Akk, UpperTriangular) + if info != 0 + return UpperTriangular(A), info + end + A[k,k] = Akk + AkkInv = inv(Akk) + for j = (k+1:n) ∩ rs + @simd for i = (1:k-1) ∩ cs ∩ colsupport(A,j) + A[i,k] -= A[i,j]*A[k,j]' + end + end + for i = (1:k-1) ∩ cs + A[i,k] *= AkkInv' + end + end + end + return UpperTriangular(A), convert(BlasInt, 0) +end + +function _reverse_chol!(A::AbstractMatrix, ::Type{LowerTriangular}) + require_one_based_indexing(A) + n = checksquare(A) + realdiag = eltype(A) <: Complex + @inbounds begin + for k = n:-1:1 + cs = colsupport(A, k) + rs = rowsupport(A, k) + Akk = realdiag ? real(A[k,k]) : A[k,k] + for i = (k+1:n) ∩ cs + Akk -= realdiag ? abs2(A[i,k]) : A[i,k]'A[i,k] + end + A[k,k] = Akk + Akk, info = _reverse_chol!(Akk, LowerTriangular) + if info != 0 + return LowerTriangular(A), info + end + A[k,k] = Akk + AkkInv = inv(copy(Akk')) + for j = (1:k-1) ∩ rs # colsupport == rowsupport + for i = (k+1:n) ∩ cs ∩ colsupport(A,j) + A[k,j] -= A[i,k]'*A[i,j] + end + A[k,j] = AkkInv*A[k,j] + end + end + end + return UpperTriangular(A), convert(BlasInt, 0) +end + +## Numbers +_reverse_chol!(x::Number, uplo) = LinearAlgebra._chol!(x, uplo) + +## for StridedMatrices, check that matrix is symmetric/Hermitian + +# cholesky!. Destructive methods for computing Cholesky factorization of real symmetric +# or Hermitian matrix +## No pivoting (default) +function reversecholesky!(A::RealHermSymComplexHerm, ::NoPivot = NoPivot(); check::Bool = true) + C, info = _reverse_chol!(A.data, A.uplo == 'U' ? UpperTriangular : LowerTriangular) + check && checkpositivedefinite(info) + return ReverseCholesky(C.data, A.uplo, info) +end + +### for AbstractMatrix, check that matrix is symmetric/Hermitian +""" + reversecholesky!(A::AbstractMatrix, NoPivot(); check = true) -> ReverseCholesky + +The same as [`reversecholesky`](@ref), but saves space by overwriting the input `A`, +instead of creating a copy. An [`InexactError`](@ref) exception is thrown if +the factorization produces a number not representable by the element type of +`A`, e.g. for integer types. +``` +""" +function reversecholesky!(A::AbstractMatrix, ::NoPivot = NoPivot(); check::Bool = true) + checksquare(A) + if !ishermitian(A) # return with info = -1 if not Hermitian + check && checkpositivedefinite(-1) + return ReverseCholesky(A, 'U', convert(BlasInt, -1)) + else + return reversecholesky!(Hermitian(A), NoPivot(); check = check) + end +end + +reversecholcopy(A) = cholcopy(A) +function reversecholcopy(A::SymTridiagonal) + T = LinearAlgebra.choltype(A) + Symmetric(Bidiagonal(AbstractVector{T}(A.dv), AbstractVector{T}(A.ev), :U)) +end + +function reversecholcopy(A::Symmetric{<:Any,<:Tridiagonal}) + T = LinearAlgebra.choltype(A) + if A.uplo == 'U' + Symmetric(Bidiagonal(AbstractVector{T}(parent(A).d), AbstractVector{T}(parent(A).du), :U)) + else + Symmetric(Bidiagonal(AbstractVector{T}(parent(A).d), AbstractVector{T}(parent(A).dl), :L), :L) + end +end + +_copyifsametype(::Type{T}, A::AbstractMatrix{T}) where T = copy(A) +_copyifsametype(_, A) = A + +function reversecholcopy(A::Symmetric{<:Any,<:Bidiagonal}) + T = LinearAlgebra.choltype(A) + B = _copyifsametype(T, AbstractMatrix{T}(parent(A))) + Symmetric{T,typeof(B)}(B, A.uplo) +end + +# reversecholesky. Non-destructive methods for computing ReverseCholesky factorization of real symmetric +# or Hermitian matrix +## No pivoting (default) +""" + reversecholesky(A, NoPivot(); check = true) -> ReverseCholesky + +Compute the ReverseCholesky factorization of a dense symmetric positive definite matrix `A` +and return a [`ReverseCholesky`](@ref) factorization. The matrix `A` can either be a [`Symmetric`](@ref) or [`Hermitian`](@ref) +[`AbstractMatrix`](@ref) or a *perfectly* symmetric or Hermitian `AbstractMatrix`. + +The triangular ReverseCholesky factor can be obtained from the factorization `F` via `F.L` and `F.U`, +where `A ≈ F.U * F.U' ≈ F.L' * F.L`. +""" +reversecholesky(A::AbstractMatrix, ::NoPivot=NoPivot(); check::Bool = true) = + reversecholesky!(reversecholcopy(A); check) + +function reversecholesky(A::AbstractMatrix{Float16}, ::NoPivot=NoPivot(); check::Bool = true) + X = reversecholesky!(reversecholcopy(A); check = check) + return ReverseCholesky{Float16}(X) +end + + +## Number +function reversecholesky(x::Number, uplo::Symbol=:U) + C, info = _reverse_chol!(x, uplo) + xf = fill(C, 1, 1) + ReverseCholesky(xf, uplo, info) +end + + +function ReverseCholesky{T}(C::ReverseCholesky) where T + Cnew = convert(AbstractMatrix{T}, C.factors) + ReverseCholesky{T, typeof(Cnew)}(Cnew, C.uplo, C.info) +end +Factorization{T}(C::ReverseCholesky{T}) where {T} = C +Factorization{T}(C::ReverseCholesky) where {T} = ReverseCholesky{T}(C) + +AbstractMatrix(C::ReverseCholesky) = C.uplo == 'U' ? C.U*C.U' : C.L'*C.L +AbstractArray(C::ReverseCholesky) = AbstractMatrix(C) +Matrix(C::ReverseCholesky) = Array(AbstractArray(C)) +Array(C::ReverseCholesky) = Matrix(C) + +copy(C::ReverseCholesky) = ReverseCholesky(copy(C.factors), C.uplo, C.info) + + +size(C::ReverseCholesky) = size(C.factors) +size(C::ReverseCholesky, d::Integer) = size(C.factors, d) + +function getproperty(C::ReverseCholesky, d::Symbol) + Cfactors = getfield(C, :factors) + Cuplo = getfield(C, :uplo) + if d === :U + return UpperTriangular(Cuplo === char_uplo(d) ? Cfactors : copy(Cfactors')) + elseif d === :L + return LowerTriangular(Cuplo === char_uplo(d) ? Cfactors : copy(Cfactors')) + elseif d === :UL + return (Cuplo === 'U' ? UpperTriangular(Cfactors) : LowerTriangular(Cfactors)) + else + return getfield(C, d) + end +end +Base.propertynames(F::ReverseCholesky, private::Bool=false) = + (:U, :L, :UL, (private ? fieldnames(typeof(F)) : ())...) + + + +issuccess(C::ReverseCholesky) = C.info == 0 + +adjoint(C::ReverseCholesky) = C + +function show(io::IO, mime::MIME{Symbol("text/plain")}, C::ReverseCholesky) + if issuccess(C) + summary(io, C); println(io) + println(io, "$(C.uplo) factor:") + show(io, mime, C.UL) + else + print(io, "Failed factorization of type $(typeof(C))") + end +end + +function ldiv!(C::ReverseCholesky, B::AbstractVecOrMat) + if C.uplo == 'L' + return ldiv!(LowerTriangular(C.factors), ldiv!(adjoint(LowerTriangular(C.factors)), B)) + else + return ldiv!(adjoint(UpperTriangular(C.factors)), ldiv!(UpperTriangular(C.factors), B)) + end +end + +function rdiv!(B::AbstractMatrix, C::ReverseCholesky) + if C.uplo == 'L' + return rdiv!(rdiv!(B, LowerTriangular(C.factors)), adjoint(LowerTriangular(C.factors))) + else + return rdiv!(rdiv!(B, adjoint(UpperTriangular(C.factors))), UpperTriangular(C.factors)) + end +end + +isposdef(C::ReverseCholesky) = C.info == 0 + +function det(C::ReverseCholesky) + dd = one(real(eltype(C))) + @inbounds for i in 1:size(C.factors,1) + dd *= real(C.factors[i,i])^2 + end + return dd +end + +function logdet(C::ReverseCholesky) + dd = zero(real(eltype(C))) + @inbounds for i in 1:size(C.factors,1) + dd += log(real(C.factors[i,i])) + end + dd + dd # instead of 2.0dd which can change the type +end + + + +logabsdet(C::ReverseCholesky) = logdet(C), one(eltype(C)) # since C is p.s.d. + + +function getproperty(C::ReverseCholesky{<:Any,<:Diagonal}, d::Symbol) + Cfactors = getfield(C, :factors) + if d in (:U, :L, :UL) + return Cfactors + else + return getfield(C, d) + end +end + +function getproperty(C::ReverseCholesky{<:Any, <:Bidiagonal}, d::Symbol) + Cfactors = getfield(C, :factors) + Cuplo = getfield(C, :uplo) + @assert Cfactors.uplo === Cuplo + if d === :U && Cuplo === 'U' + return Cfactors + elseif d === :L && Cuplo === 'U' + return Cfactors' + elseif d === :U && Cuplo === 'L' + return Cfactors' + elseif d === :L && Cuplo === 'L' + return Cfactors + elseif d === :UL + return Cfactors + else + return getfield(C, d) + end +end + +inv(C::ReverseCholesky{<:Any,<:Diagonal}) = Diagonal(map(inv∘abs2, C.factors.diag)) \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index b092d32..77d1c66 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -439,33 +439,5 @@ end end include("test_rq.jl") - -# choleskyinv.jl -@testset "choleskyinv" begin - etol = 1e-9 - n = 20 - # real matrices - A = randn(n, n) - P = A*A' - Pi = inv(P) - C = choleskyinv!(copy(Matrix(P))) - @test(norm(C.c.L*C.c.U-P)/√n < etol) - @test(norm(C.ci.U*C.ci.L-Pi)/√n < etol) - C = choleskyinv(P) - @test(norm(C.c.L*C.c.U-P)/√n < etol) - @test(norm(C.ci.U*C.ci.L-Pi)/√n < etol) - - # repeat the test for complex matrices - A=randn(ComplexF64, n, n) - P=A*A' - Pi=inv(P) - C=choleskyinv!(copy(Matrix(P))) - @test(norm(C.c.L*C.c.U-P)/√n < etol) - @test(norm(C.ci.U*C.ci.L-Pi)/√n < 4etol) - C = choleskyinv(P) - @test(norm(C.c.L*C.c.U-P)/√n < etol) - @test(norm(C.ci.U*C.ci.L-Pi)/√n < 4etol) -end - - include("test_polar.jl") +include("test_reversecholesky.jl") \ No newline at end of file diff --git a/test/test_reversecholesky.jl b/test/test_reversecholesky.jl new file mode 100644 index 0000000..9be9c7b --- /dev/null +++ b/test/test_reversecholesky.jl @@ -0,0 +1,415 @@ +# This file is modifed from a part of Julia. License is MIT: https://julialang.org/license + +using Test, LinearAlgebra, Random, MatrixFactorizations +using LinearAlgebra: BlasComplex, BlasFloat, BlasReal, QRPivoted, + PosDefException, RankDeficientException, chkfullrank + +function unary_ops_tests(a, ca, tol; n=size(a, 1)) + @test inv(ca)*a ≈ Matrix(I, n, n) + @test a*inv(ca) ≈ Matrix(I, n, n) + @test abs((det(ca) - det(a))/det(ca)) <= tol # Ad hoc, but statistically verified, revisit + @test logdet(ca) ≈ logdet(a) + @test logdet(ca) ≈ log(det(ca)) # logdet is less likely to overflow + logabsdet_ca = logabsdet(ca) + logabsdet_a = logabsdet(a) + @test logabsdet_ca[1] ≈ logabsdet_a[1] + @test logabsdet_ca[2] ≈ logabsdet_a[2] + @test isposdef(ca) + @test_throws ErrorException ca.Z + @test size(ca) == size(a) + @test Array(copy(ca)) ≈ a +end + +function factor_recreation_tests(a_U, a_L) + c_U = reversecholesky(a_U) + c_L = reversecholesky(a_L) + cl = c_L.U + ls = c_L.L + @test Array(c_U) ≈ Array(c_L) ≈ a_U + @test ls'*ls ≈ a_U + @test triu(c_U.factors) ≈ c_U.U + @test tril(c_L.factors) ≈ c_L.L + @test istriu(cl) + @test cl*cl' ≈ a_U + @test cl*cl' ≈ a_L +end + +@testset "ReverseCholesky" begin + @testset "core functionality" begin + n = 10 + + # Split n into 2 parts for tests needing two matrices + n1 = div(n, 2) + n2 = 2*n1 + + Random.seed!(12344) + + areal = randn(n,n)/2 + aimg = randn(n,n)/2 + a2real = randn(n,n)/2 + a2img = randn(n,n)/2 + breal = randn(n,2)/2 + bimg = randn(n,2)/2 + + @testset "basic sanity check" begin + A = [4 2; 2 4] + U = reversecholesky(A).U + @test U*U' ≈ A + @test cholesky(A).L[end:-1:1,end:-1:1] ≈ U + + A = Symmetric(areal) + 10I + U = reversecholesky(A).U + @test U*U' ≈ A + + A = Hermitian(areal + im*aimg) + 10I + U = reversecholesky(A).U + @test U*U' ≈ A + + A = Symmetric(areal, :L) + 10I + U = reversecholesky(A).U + @test U*U' ≈ A + + A = Hermitian(areal + im*aimg, :L) + 10I + U = reversecholesky(A).U + @test U*U' ≈ A + end + + for eltya in (Float32, Float64, ComplexF32, ComplexF64, BigFloat, Int) + a = eltya == Int ? rand(1:7, n, n) : convert(Matrix{eltya}, eltya <: Complex ? complex.(areal, aimg) : areal) + a2 = eltya == Int ? rand(1:7, n, n) : convert(Matrix{eltya}, eltya <: Complex ? complex.(a2real, a2img) : a2real) + + ε = εa = eps(abs(float(one(eltya)))) + + # Test of symmetric pos. def. strided matrix + apd = a'*a + @inferred reversecholesky(apd) + capd = reversecholesky(apd) + r = capd.U + κ = cond(apd, 1) #condition number + + unary_ops_tests(apd, capd, ε*κ*n) + if eltya != Int + @test Factorization{eltya}(capd) === capd + if eltya <: Real + @test Array(Factorization{complex(eltya)}(capd)) ≈ Array(factorize(complex(apd))) + @test eltype(Factorization{complex(eltya)}(capd)) == complex(eltya) + end + end + @testset "throw for non-square input" begin + A = rand(eltya, 2, 3) + @test_throws DimensionMismatch reversecholesky(A) + @test_throws DimensionMismatch reversecholesky!(A) + end + + #Test error bound on reconstruction of matrix: LAWNS 14, Lemma 2.1 + + #these tests were failing on 64-bit linux when inside the inner loop + #for eltya = ComplexF32 and eltyb = Int. The E[i,j] had NaN32 elements + #but only with Random.seed!(1234321) set before the loops. + E = abs.(apd - r*r') + for i=1:n, j=1:n + @test E[i,j] <= (n+1)ε/(1-(n+1)ε)*real(sqrt(apd[i,i]*apd[j,j])) + end + E = abs.(apd - Matrix(capd)) + for i=1:n, j=1:n + @test E[i,j] <= (n+1)ε/(1-(n+1)ε)*real(sqrt(apd[i,i]*apd[j,j])) + end + @test LinearAlgebra.issuccess(capd) + @inferred(logdet(capd)) + + apos = apd[1,1] + @test all(x -> x ≈ √apos, cholesky(apos).factors) + + # Test cholesky with Symmetric/Hermitian upper/lower + apds = Symmetric(apd) + apdsL = Symmetric(apd, :L) + apdh = Hermitian(apd) + apdhL = Hermitian(apd, :L) + if eltya <: Real + capds = reversecholesky(apds) + unary_ops_tests(apds, capds, ε*κ*n) + if eltya <: BlasReal + capds = cholesky!(copy(apds)) + unary_ops_tests(apds, capds, ε*κ*n) + end + ulstring = sprint((t, s) -> show(t, "text/plain", s), capds.UL) + @test sprint((t, s) -> show(t, "text/plain", s), capds) == "$(typeof(capds))\nU factor:\n$ulstring" + else + capdh = reversecholesky(apdh) + unary_ops_tests(apdh, capdh, ε*κ*n) + capdh = reversecholesky!(copy(apdh)) + unary_ops_tests(apdh, capdh, ε*κ*n) + capdh = reversecholesky!(copy(apd)) + unary_ops_tests(apd, capdh, ε*κ*n) + ulstring = sprint((t, s) -> show(t, "text/plain", s), capdh.UL) + @test sprint((t, s) -> show(t, "text/plain", s), capdh) == "$(typeof(capdh))\nU factor:\n$ulstring" + end + + # test reversecholesky of 2x2 Strang matrix + S = SymTridiagonal{eltya}([2, 2], [-1]) + for uplo in (:U, :L) + @test Matrix(@inferred reversecholesky(Hermitian(S, uplo))) ≈ S + if eltya <: Real + @test Matrix(@inferred reversecholesky(Symmetric(S, uplo))) ≈ S + end + end + @test Matrix(reversecholesky(S).U) ≈ [sqrt(eltya(3)) -1; 0 2] / sqrt(eltya(2)) + @test Matrix(reversecholesky(S)) ≈ S + + # test extraction of factor and re-creating original matrix + if eltya <: Real + factor_recreation_tests(apds, apdsL) + else + factor_recreation_tests(apdh, apdhL) + end + + + + for eltyb in (Float32, Float64, ComplexF32, ComplexF64, Int) + b = eltyb == Int ? rand(1:5, n, 2) : convert(Matrix{eltyb}, eltyb <: Complex ? complex.(breal, bimg) : breal) + εb = eps(abs(float(one(eltyb)))) + ε = max(εa,εb) + + for b in (b, view(b, 1:n, 1)) # Array and SubArray + + # Test error bound on linear solver: LAWNS 14, Theorem 2.1 + # This is a surprisingly loose bound + x = capd\b + @test norm(x-apd\b,1)/norm(x,1) <= (3n^2 + n + n^3*ε)*ε/(1-(n+1)*ε)*κ + @test norm(apd*x-b,1)/norm(b,1) <= (3n^2 + n + n^3*ε)*ε/(1-(n+1)*ε)*κ + + @test norm(a*(capd\(a'*b)) - b,1)/norm(b,1) <= ε*κ*n # Ad hoc, revisit + + if eltya != BigFloat && eltyb != BigFloat + lapd = reversecholesky(apdhL) + @test norm(apd * (lapd\b) - b)/norm(b) <= ε*κ*n + @test norm(apd * (lapd\b[1:n]) - b[1:n])/norm(b[1:n]) <= ε*κ*n + end + end + end + + for eltyb in (Float64, ComplexF64) + Breal = convert(Matrix{BigFloat}, randn(n,n)/2) + Bimg = convert(Matrix{BigFloat}, randn(n,n)/2) + B = (eltya <: Complex || eltyb <: Complex) ? complex.(Breal, Bimg) : Breal + εb = eps(abs(float(one(eltyb)))) + ε = max(εa,εb) + + for B in (B, view(B, 1:n, 1:n)) # Array and SubArray + + # Test error bound on linear solver: LAWNS 14, Theorem 2.1 + # This is a surprisingly loose bound + BB = copy(B) + ldiv!(capd, BB) + @test norm(apd \ B - BB, 1) / norm(BB, 1) <= (3n^2 + n + n^3*ε)*ε/(1-(n+1)*ε)*κ + @test norm(apd * BB - B, 1) / norm(B, 1) <= (3n^2 + n + n^3*ε)*ε/(1-(n+1)*ε)*κ + end + end + + @testset "solve with generic Cholesky" begin + Breal = convert(Matrix{BigFloat}, randn(n,n)/2) + Bimg = convert(Matrix{BigFloat}, randn(n,n)/2) + B = eltya <: Complex ? complex.(Breal, Bimg) : Breal + εb = eps(abs(float(one(eltype(B))))) + ε = max(εa,εb) + + for B in (B, view(B, 1:n, 1:n)) # Array and SubArray + + # Test error bound on linear solver: LAWNS 14, Theorem 2.1 + # This is a surprisingly loose bound + cpapd = reversecholesky(eltya <: Complex ? apdh : apds) + BB = copy(B) + rdiv!(BB, cpapd) + @test norm(B / apd - BB, 1) / norm(BB, 1) <= (3n^2 + n + n^3*ε)*ε/(1-(n+1)*ε)*κ + @test norm(BB * apd - B, 1) / norm(B, 1) <= (3n^2 + n + n^3*ε)*ε/(1-(n+1)*ε)*κ + cpapd = reversecholesky(eltya <: Complex ? apdhL : apdsL) + BB = copy(B) + rdiv!(BB, cpapd) + @test norm(B / apd - BB, 1) / norm(BB, 1) <= (3n^2 + n + n^3*ε)*ε/(1-(n+1)*ε)*κ + @test norm(BB * apd - B, 1) / norm(B, 1) <= (3n^2 + n + n^3*ε)*ε/(1-(n+1)*ε)*κ + end + end + end + end + + @testset "behavior for non-positive definite matrices" for T in (Float64, ComplexF64) + A = T[1 2; 2 1] + B = T[1 2; 0 1] + # check = (true|false) + for M in (A, Hermitian(A), B) + @test_throws PosDefException reversecholesky(M) + @test_throws PosDefException reversecholesky!(copy(M)) + @test_throws PosDefException reversecholesky(M; check = true) + @test_throws PosDefException reversecholesky!(copy(M); check = true) + @test !LinearAlgebra.issuccess(reversecholesky(M; check = false)) + @test !LinearAlgebra.issuccess(reversecholesky!(copy(M); check = false)) + end + str = sprint((io, x) -> show(io, "text/plain", x), reversecholesky(A; check = false)) + end + + @testset "Cholesky factor of Matrix with non-commutative elements, here 2x2-matrices" begin + X = Matrix{Float64}[0.1*rand(2,2) for i in 1:3, j = 1:3] + L = Matrix(MatrixFactorizations._reverse_chol!(X*X', LowerTriangular)[1]) + U = Matrix(MatrixFactorizations._reverse_chol!(X*X', UpperTriangular)[1]) + XX = Matrix(X*X') + + @test_broken sum(sum(norm, L'*L - XX)) < eps() + @test_broken sum(sum(norm, U*U' - XX)) < eps() + end + + @testset "Non-strided Cholesky solves" begin + B = randn(5, 5) + v = rand(5) + @test reversecholesky(Diagonal(v)) \ B ≈ Diagonal(v) \ B + @test B / reversecholesky(Diagonal(v)) ≈ B / Diagonal(v) + @test inv(reversecholesky(Diagonal(v)))::Diagonal ≈ Diagonal(1 ./ v) + end + + + @testset "cholesky Diagonal" begin + # real + d = abs.(randn(3)) .+ 0.1 + D = Diagonal(d) + CD = reversecholesky(D) + CM = reversecholesky(Matrix(D)) + @test CD isa ReverseCholesky{Float64} + @test CD.U ≈ Diagonal(.√d) ≈ CM.U + @test D ≈ CD.L * CD.U + @test CD.info == 0 + + F = reversecholesky(Hermitian(I(3))) + @test F isa ReverseCholesky{Float64,<:Diagonal} + @test Matrix(F) ≈ I(3) + + # real, failing + @test_throws PosDefException reversecholesky(Diagonal([1.0, -2.0])) + Dnpd = reversecholesky(Diagonal([1.0, -2.0]); check = false) + @test Dnpd.info == 1 + + # complex + D = complex(D) + CD = reversecholesky(Hermitian(D)) + CM = reversecholesky(Matrix(Hermitian(D))) + @test CD isa ReverseCholesky{ComplexF64,<:Diagonal} + @test CD.U ≈ Diagonal(.√d) ≈ CM.U + @test D ≈ CD.L * CD.U + @test CD.info == 0 + + # complex, failing + D[2, 2] = 0.0 + 0im + @test reversecholesky(D).U ≈ sqrt.(D) + + + # InexactError for Int + @test_throws InexactError reversecholesky!(Diagonal([2, 1])) + end + + @testset "Cholesky for AbstractMatrix" begin + S = SymTridiagonal(fill(2.0, 4), ones(3)) + C = reversecholesky(S) + @test C.U * C.L ≈ S + end + + @testset "constructor with non-BlasInt arguments" begin + x = rand(5,5) + chol = reversecholesky(x'x) + + factors, uplo, info = chol.factors, chol.uplo, chol.info + + @test ReverseCholesky(factors, uplo, Int32(info)) == chol + @test ReverseCholesky(factors, uplo, Int64(info)) == chol + end + + @testset "issue #37356, diagonal elements of hermitian generic matrix" begin + B = Hermitian(hcat([one(BigFloat) + im])) + @test Matrix(reversecholesky(B)) ≈ B + C = Hermitian(hcat([one(BigFloat) + im]), :L) + @test Matrix(reversecholesky(C)) ≈ C + end + + @testset "constructing a ReverseCholesky factor from a triangular matrix" begin + A = [1.0 2.0; 3.0 4.0] + let + U = UpperTriangular(A) + C = ReverseCholesky(U) + @test C isa ReverseCholesky{Float64} + @test C.U == U + @test C.L == U' + end + let + L = LowerTriangular(A) + C = ReverseCholesky(L) + @test C isa ReverseCholesky{Float64} + @test C.L == L + @test C.U == L' + end + end + + @testset "adjoint of ReverseCholesky" begin + A = randn(5, 5) + A = A'A + F = reversecholesky(A) + b = ones(size(A, 1)) + @test F\b == F'\b + end + + @testset "Float16" begin + A = Float16[4. 12. -16.; 12. 37. -43.; -16. -43. 98.] + B = reversecholesky(A) + B32 = reversecholesky(Float32.(A)) + @test B isa ReverseCholesky{Float16, Matrix{Float16}} + @test B.U isa UpperTriangular{Float16, Matrix{Float16}} + @test B.L isa LowerTriangular{Float16, Matrix{Float16}} + @test B.UL isa UpperTriangular{Float16, Matrix{Float16}} + @test B.U ≈ B32.U + @test B.L ≈ B32.L + @test B.UL ≈ B32.UL + @test Matrix(B) ≈ A + end + + @testset "large Sparse" begin + n = 1_000_000 + A = SymTridiagonal(fill(4,n), fill(1,n-1)) + R = reversecholesky(A) + @test R isa ReverseCholesky{Float64, <:Bidiagonal} + @test R.U isa Bidiagonal + @test R.L isa Bidiagonal + @test R.UL isa Bidiagonal + # Bidiagonal multiplication not supported + @test R.U*(R.U' * [1; zeros(n-1)]) ≈ A[:,1] + @test R.L'*(R.L * [1; zeros(n-1)]) ≈ A[:,1] + + A = Tridiagonal(fill(1.0,n-1), fill(4.0,n), fill(1/2,n-1)) + R = reversecholesky(Symmetric(A)) + @test R.U*(R.U' * [1; zeros(n-1)]) ≈ A[1,:] + @test R.L'*(R.L * [1; zeros(n-1)]) ≈ A[1,:] + R = reversecholesky(Symmetric(A,:L)) + @test R.U*(R.U' * [1; zeros(n-1)]) ≈ A[:,1] + @test R.L'*(R.L * [1; zeros(n-1)]) ≈ A[:,1] + + A = Bidiagonal(fill(4.0,n), fill(1.0,n-1), :U) + R = reversecholesky(Symmetric(A)) + @test R.U*(R.U' * [1; zeros(n-1)]) ≈ A[1,:] + @test R.L'*(R.L * [1; zeros(n-1)]) ≈ A[1,:] + + A = Bidiagonal(fill(4.0,n), fill(1.0,n-1), :L) + R = reversecholesky(Symmetric(A,:L)) + @test R.U * (R.U' * [1; zeros(n-1)]) ≈ A[:,1] + @test R.L'*(R.L * [1; zeros(n-1)]) ≈ A[:,1] + + A = Bidiagonal(fill(4,n), fill(1,n-1), :L) + R = reversecholesky(Symmetric(A,:L)) + @test R.U * (R.U' * [1; zeros(n-1)]) ≈ A[:,1] + @test R.L'*(R.L * [1; zeros(n-1)]) ≈ A[:,1] + end + + @testset "coverage" begin + A = [4 2; 2 4] + R = reversecholesky(A) + @test ReverseCholesky(R.factors, R.uplo, R.info) == R + U,L = R + @test U*U' ≈ A + @test reversecholesky(5).U[1,1] ≈ sqrt(5) + @test propertynames(R) == (:U, :L, :UL) + end +end \ No newline at end of file