use testsets in cholesky test (#22138)

(cherry picked from commit b11bc92)
JuliaLang · Sep 12, 2017 · 96582db · 96582db
1 parent 7d088ba
commit 96582db
Showing 1 changed file with 168 additions and 164 deletions.
diff --git a/test/linalg/cholesky.jl b/test/linalg/cholesky.jl
@@ -1,165 +1,157 @@
 # This file is a part of Julia. License is MIT: https://julialang.org/license
 
-debug = false
-
 using Base.Test
 
 using Base.LinAlg: BlasComplex, BlasFloat, BlasReal, QRPivoted
 
-n = 10
+@testset "core functionality" begin
+    n = 10
 
-# Split n into 2 parts for tests needing two matrices
-n1 = div(n, 2)
-n2 = 2*n1
+    # Split n into 2 parts for tests needing two matrices
+    n1 = div(n, 2)
+    n2 = 2*n1
 
-srand(1234321)
+    srand(1234321)
 
-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
+    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
 
-for eltya in (Float32, Float64, Complex64, Complex128, 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)
-    apd  = a'*a                  # symmetric positive-definite
+    for eltya in (Float32, Float64, Complex64, Complex128, 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)
+        apd  = a'*a                  # symmetric positive-definite
 
-    apds  = Symmetric(apd)
-    apdsL = Symmetric(apd, :L)
-    apdh  = Hermitian(apd)
-    apdhL = Hermitian(apd, :L)
-    ε = εa = eps(abs(float(one(eltya))))
+        apds  = Symmetric(apd)
+        apdsL = Symmetric(apd, :L)
+        apdh  = Hermitian(apd)
+        apdhL = Hermitian(apd, :L)
+        ε = εa = eps(abs(float(one(eltya))))
 
-    @inferred cholfact(apd)
-    @inferred chol(apd)
-    capd  = factorize(apd)
-    r     = capd[:U]
-    κ     = cond(apd, 1) #condition number
+        @inferred cholfact(apd)
+        @inferred chol(apd)
+        capd  = factorize(apd)
+        r     = capd[:U]
+        κ     = cond(apd, 1) #condition number
 
-    #getindex
-    @test_throws KeyError capd[:Z]
+        #getindex
+        @test_throws KeyError capd[:Z]
 
-    #Test error bound on reconstruction of matrix: LAWNS 14, Lemma 2.1
+        #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 = Complex64 and eltyb = Int. The E[i,j] had NaN32 elements
-    #but only with srand(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 - full(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 apd*inv(capd) ≈ eye(n)
-    @test abs((det(capd) - det(apd))/det(capd)) <= ε*κ*n # Ad hoc, but statistically verified, revisit
-    @test @inferred(logdet(capd)) ≈ log(det(capd)) # logdet is less likely to overflow
-
-    apos = apd[1,1]            # test chol(x::Number), needs x>0
-    @test all(x -> x ≈ √apos, cholfact(apos).factors)
-    @test_throws ArgumentError chol(-one(eltya))
-
-    if eltya <: Real
-        capds = cholfact(apds)
-        @test inv(capds)*apds ≈ eye(n)
-        @test abs((det(capds) - det(apd))/det(capds)) <= ε*κ*n
-        if eltya <: BlasReal
-            capds = cholfact!(copy(apds))
+        #these tests were failing on 64-bit linux when inside the inner loop
+        #for eltya = Complex64 and eltyb = Int. The E[i,j] had NaN32 elements
+        #but only with srand(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 - full(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 apd*inv(capd) ≈ eye(n)
+        @test abs((det(capd) - det(apd))/det(capd)) <= ε*κ*n # Ad hoc, but statistically verified, revisit
+        @test @inferred(logdet(capd)) ≈ log(det(capd)) # logdet is less likely to overflow
+
+        apos = apd[1,1]            # test chol(x::Number), needs x>0
+        @test all(x -> x ≈ √apos, cholfact(apos).factors)
+        @test_throws PosDefException chol(-one(eltya))
+
+        if eltya <: Real
+            capds = cholfact(apds)
             @test inv(capds)*apds ≈ eye(n)
             @test abs((det(capds) - det(apd))/det(capds)) <= ε*κ*n
+            if eltya <: BlasReal
+                capds = cholfact!(copy(apds))
+                @test inv(capds)*apds ≈ eye(n)
+                @test abs((det(capds) - det(apd))/det(capds)) <= ε*κ*n
+            end
+            ulstring = sprint(show,capds[:UL])
+            @test sprint(show,capds) == "$(typeof(capds)) with factor:\n$ulstring"
+        else
+            capdh = cholfact(apdh)
+            @test inv(capdh)*apdh ≈ eye(n)
+            @test abs((det(capdh) - det(apd))/det(capdh)) <= ε*κ*n
+            capdh = cholfact!(copy(apdh))
+            @test inv(capdh)*apdh ≈ eye(n)
+            @test abs((det(capdh) - det(apd))/det(capdh)) <= ε*κ*n
+            capdh = cholfact!(copy(apd))
+            @test inv(capdh)*apdh ≈ eye(n)
+            @test abs((det(capdh) - det(apd))/det(capdh)) <= ε*κ*n
+            capdh = cholfact!(copy(apd), :L)
+            @test inv(capdh)*apdh ≈ eye(n)
+            @test abs((det(capdh) - det(apd))/det(capdh)) <= ε*κ*n
+            ulstring = sprint(show,capdh[:UL])
+            @test sprint(show,capdh) == "$(typeof(capdh)) with factor:\n$ulstring"
         end
-        ulstring = sprint(show,capds[:UL])
-        @test sprint(show,capds) == "$(typeof(capds)) with factor:\n$ulstring"
-    else
-        capdh = cholfact(apdh)
-        @test inv(capdh)*apdh ≈ eye(n)
-        @test abs((det(capdh) - det(apd))/det(capdh)) <= ε*κ*n
-        capdh = cholfact!(copy(apdh))
-        @test inv(capdh)*apdh ≈ eye(n)
-        @test abs((det(capdh) - det(apd))/det(capdh)) <= ε*κ*n
-        capdh = cholfact!(copy(apd))
-        @test inv(capdh)*apdh ≈ eye(n)
-        @test abs((det(capdh) - det(apd))/det(capdh)) <= ε*κ*n
-        capdh = cholfact!(copy(apd), :L)
-        @test inv(capdh)*apdh ≈ eye(n)
-        @test abs((det(capdh) - det(apd))/det(capdh)) <= ε*κ*n
-        ulstring = sprint(show,capdh[:UL])
-        @test sprint(show,capdh) == "$(typeof(capdh)) with factor:\n$ulstring"
-    end
 
-    # test chol of 2x2 Strang matrix
-    S = convert(AbstractMatrix{eltya},full(SymTridiagonal([2,2],[-1])))
-    U = Bidiagonal([2,sqrt(eltya(3))],[-1],true) / sqrt(eltya(2))
-    @test full(chol(S)) ≈ full(U)
-
-    #lower Cholesky factor
-    lapd = cholfact(apd, :L)
-    @test full(lapd) ≈ apd
-    l = lapd[:L]
-    @test l*l' ≈ apd
-    @test triu(capd.factors) ≈ lapd[:U]
-    @test tril(lapd.factors) ≈ capd[:L]
-    if eltya <: Real
-        capds = cholfact(apds)
-        lapds = cholfact(apdsL)
-        cl    = chol(apdsL)
-        ls = lapds[:L]
-        @test ls*ls' ≈ apd
-        @test triu(capds.factors) ≈ lapds[:U]
-        @test tril(lapds.factors) ≈ capds[:L]
-        @test istriu(cl)
-        @test cl'cl ≈ apds
-        @test cl'cl ≈ apdsL
-    else
-        capdh = cholfact(apdh)
-        lapdh = cholfact(apdhL)
-        cl    = chol(apdhL)
-        ls = lapdh[:L]
-        @test ls*ls' ≈ apd
-        @test triu(capdh.factors) ≈ lapdh[:U]
-        @test tril(lapdh.factors) ≈ capdh[:L]
-        @test istriu(cl)
-        @test cl'cl ≈ apdh
-        @test cl'cl ≈ apdhL
-    end
+        # test chol of 2x2 Strang matrix
+        S = convert(AbstractMatrix{eltya},full(SymTridiagonal([2,2],[-1])))
+        U = Bidiagonal([2,sqrt(eltya(3))],[-1],true) / sqrt(eltya(2))
+        @test full(chol(S)) ≈ full(U)
+
+        #lower Cholesky factor
+        lapd = cholfact(apd, :L)
+        @test full(lapd) ≈ apd
+        l = lapd[:L]
+        @test l*l' ≈ apd
+        @test triu(capd.factors) ≈ lapd[:U]
+        @test tril(lapd.factors) ≈ capd[:L]
+        if eltya <: Real
+            capds = cholfact(apds)
+            lapds = cholfact(apdsL)
+            cl    = chol(apdsL)
+            ls = lapds[:L]
+            @test ls*ls' ≈ apd
+            @test triu(capds.factors) ≈ lapds[:U]
+            @test tril(lapds.factors) ≈ capds[:L]
+            @test istriu(cl)
+            @test cl'cl ≈ apds
+            @test cl'cl ≈ apdsL
+        else
+            capdh = cholfact(apdh)
+            lapdh = cholfact(apdhL)
+            cl    = chol(apdhL)
+            ls = lapdh[:L]
+            @test ls*ls' ≈ apd
+            @test triu(capdh.factors) ≈ lapdh[:U]
+            @test tril(lapdh.factors) ≈ capdh[:L]
+            @test istriu(cl)
+            @test cl'cl ≈ apdh
+            @test cl'cl ≈ apdhL
+        end
 
-    #pivoted upper Cholesky
-    if eltya != BigFloat
-        cz = cholfact(zeros(eltya,n,n), :U, Val{true})
-        @test_throws Base.LinAlg.RankDeficientException Base.LinAlg.chkfullrank(cz)
-        cpapd = cholfact(apd, :U, Val{true})
-        @test rank(cpapd) == n
-        @test all(diff(diag(real(cpapd.factors))).<=0.) # diagonal should be non-increasing
-        if isreal(apd)
-            @test apd*inv(cpapd) ≈ eye(n)
+        #pivoted upper Cholesky
+        if eltya != BigFloat
+            cz = cholfact(zeros(eltya,n,n), :U, Val{true})
+            @test_throws Base.LinAlg.RankDeficientException Base.LinAlg.chkfullrank(cz)
+            cpapd = cholfact(apd, :U, Val{true})
+            @test rank(cpapd) == n
+            @test all(diff(diag(real(cpapd.factors))).<=0.) # diagonal should be non-increasing
+            if isreal(apd)
+                @test apd*inv(cpapd) ≈ eye(n)
+            end
+            @test full(cpapd) ≈ apd
+            #getindex
+            @test_throws KeyError cpapd[:Z]
+
+            @test size(cpapd) == size(apd)
+            @test full(copy(cpapd)) ≈ apd
+            @test det(cpapd) ≈ det(apd)
+            @test logdet(cpapd) ≈ logdet(apd)
+            @test cpapd[:P]*cpapd[:L]*cpapd[:U]*cpapd[:P]' ≈ apd
         end
-        @test full(cpapd) ≈ apd
-        #getindex
-        @test_throws KeyError cpapd[:Z]
 
-        @test size(cpapd) == size(apd)
-        @test full(copy(cpapd)) ≈ apd
-        @test det(cpapd) ≈ det(apd)
-        @test logdet(cpapd) ≈ logdet(apd)
-        @test cpapd[:P]*cpapd[:L]*cpapd[:U]*cpapd[:P]' ≈ apd
-    end
+        for eltyb in (Float32, Float64, Complex64, Complex128, 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 eltyb in (Float32, Float64, Complex64, Complex128, 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)
-
-debug && println("\ntype of a: ", eltya, " type of b: ", eltyb, "\n")
-        let Bs = b
-            for atype in ("Array", "SubArray")
-                if atype == "Array"
-                    b = Bs
-                else
-                    b = view(Bs, 1:n, 1)
-                end
+            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
@@ -175,7 +167,6 @@ debug && println("\ntype of a: ", eltya, " type of b: ", eltyb, "\n")
                 end
                 @test_throws DimensionMismatch lapd\RowVector(ones(n))
 
-debug && println("pivoted Cholesky decomposition")
                 if eltya != BigFloat && eltyb != BigFloat # Note! Need to implement pivoted Cholesky decomposition in julia
 
                     @test norm(apd * (cpapd\b) - b)/norm(b) <= ε*κ*n # Ad hoc, revisit
@@ -192,9 +183,7 @@ debug && println("pivoted Cholesky decomposition")
     end
 end
 
-begin
-    # Cholesky factor of Matrix with non-commutative elements, here 2x2-matrices
-
+@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 = full(Base.LinAlg._chol!(X*X', LowerTriangular))
     U = full(Base.LinAlg._chol!(X*X', UpperTriangular))
@@ -205,19 +194,22 @@ begin
 end
 
 # Test generic cholfact!
-for elty in (Float32, Float64, Complex{Float32}, Complex{Float64})
-    if elty <: Complex
-        A = complex.(randn(5,5), randn(5,5))
-    else
-        A = randn(5,5)
+@testset "generic cholfact!" begin
+    for elty in (Float32, Float64, Complex{Float32}, Complex{Float64})
+        if elty <: Complex
+            A = complex.(randn(5,5), randn(5,5))
+        else
+            A = randn(5,5)
+        end
+        A = convert(Matrix{elty}, A'A)
+        @test full(cholfact(A)[:L]) ≈ full(invoke(Base.LinAlg._chol!, Tuple{AbstractMatrix, Type{LowerTriangular}}, copy(A), LowerTriangular)[1])
+        @test full(cholfact(A)[:U]) ≈ full(invoke(Base.LinAlg._chol!, Tuple{AbstractMatrix, Type{UpperTriangular}}, copy(A), UpperTriangular)[1])
     end
-    A = convert(Matrix{elty}, A'A)
-    @test full(cholfact(A)[:L]) ≈ full(invoke(Base.LinAlg._chol!, Tuple{AbstractMatrix, Type{LowerTriangular}}, copy(A), LowerTriangular))
-    @test full(cholfact(A)[:U]) ≈ full(invoke(Base.LinAlg._chol!, Tuple{AbstractMatrix, Type{UpperTriangular}}, copy(A), UpperTriangular))
 end
 
-# Test up- and downdates
-let A = complex.(randn(10,5), randn(10, 5)), v = complex.(randn(5), randn(5))
+@testset "cholesky up- and downdates" begin
+    A = complex.(randn(10,5), randn(10, 5))
+    v = complex.(randn(5), randn(5))
     for uplo in (:U, :L)
         AcA = A'A
         BcB = AcA + v*v'
@@ -231,8 +223,8 @@ let A = complex.(randn(10,5), randn(10, 5)), v = complex.(randn(5), randn(5))
     end
 end
 
-# issue #13243, unexpected nans in complex cholfact
-let apd = [5.8525753f0 + 0.0f0im -0.79540455f0 + 0.7066077f0im 0.98274714f0 + 1.3824869f0im 2.619998f0 + 1.8532984f0im -1.8306153f0 - 1.2336911f0im 0.32275113f0 + 0.015575029f0im 2.1968813f0 + 1.0640624f0im 0.27894387f0 + 0.97911835f0im 3.0476584f0 + 0.18548489f0im 0.3842994f0 + 0.7050991f0im
+@testset "issue #13243, unexpected nans in complex cholfact" begin
+    apd = [5.8525753f0 + 0.0f0im -0.79540455f0 + 0.7066077f0im 0.98274714f0 + 1.3824869f0im 2.619998f0 + 1.8532984f0im -1.8306153f0 - 1.2336911f0im 0.32275113f0 + 0.015575029f0im 2.1968813f0 + 1.0640624f0im 0.27894387f0 + 0.97911835f0im 3.0476584f0 + 0.18548489f0im 0.3842994f0 + 0.7050991f0im
         -0.79540455f0 - 0.7066077f0im 8.313246f0 + 0.0f0im -1.8076122f0 - 0.8882447f0im 0.47806996f0 + 0.48494184f0im 0.5096429f0 - 0.5395974f0im -0.7285097f0 - 0.10360408f0im -1.1760061f0 - 2.7146957f0im -0.4271084f0 + 0.042899966f0im -1.7228563f0 + 2.8335886f0im 1.8942566f0 + 0.6389735f0im
         0.98274714f0 - 1.3824869f0im -1.8076122f0 + 0.8882447f0im 9.367975f0 + 0.0f0im -0.1838578f0 + 0.6468568f0im -1.8338387f0 + 0.7064959f0im 0.041852742f0 - 0.6556877f0im 2.5673025f0 + 1.9732997f0im -1.1148382f0 - 0.15693812f0im 2.4704504f0 - 1.0389464f0im 1.0858271f0 - 1.298006f0im
         2.619998f0 - 1.8532984f0im 0.47806996f0 - 0.48494184f0im -0.1838578f0 - 0.6468568f0im 3.1117508f0 + 0.0f0im -1.956626f0 + 0.22825956f0im 0.07081801f0 - 0.31801307f0im 0.3698375f0 - 0.5400855f0im 0.80686307f0 + 1.5315914f0im 1.5649154f0 - 1.6229297f0im -0.112077385f0 + 1.2014246f0im
@@ -262,13 +254,25 @@ let apd = [5.8525753f0 + 0.0f0im -0.79540455f0 + 0.7066077f0im 0.98274714f0 + 1.
     end
 end
 
-# Fail if non-Hermitian
-@test_throws ArgumentError cholfact(randn(5,5))
-@test_throws ArgumentError cholfact(complex.(randn(5,5), randn(5,5)))
-@test_throws ArgumentError Base.LinAlg.chol!(randn(5,5))
-@test_throws ArgumentError Base.LinAlg.cholfact!(randn(5,5),:U,Val{false})
-@test_throws ArgumentError Base.LinAlg.cholfact!(randn(5,5),:U,Val{true})
-@test_throws ArgumentError cholfact(randn(5,5),:U,Val{false})
+@testset "throw if non-Hermitian" begin
+    @test_throws ArgumentError cholfact(randn(5,5))
+    @test_throws ArgumentError cholfact(complex.(randn(5,5), randn(5,5)))
+    @test_throws ArgumentError Base.LinAlg.chol!(randn(5,5))
+    @test_throws ArgumentError Base.LinAlg.cholfact!(randn(5,5),:U,Val{false})
+    @test_throws ArgumentError Base.LinAlg.cholfact!(randn(5,5),:U,Val{true})
+    @test_throws ArgumentError cholfact(randn(5,5),:U,Val{false})
+end
 
-# Fail for non-BLAS element types
-@test_throws ArgumentError cholfact!(Hermitian(rand(Float16, 5,5)), Val{true})
+@testset "fail for non-BLAS element types" begin
+    @test_throws ArgumentError cholfact!(Hermitian(rand(Float16, 5,5)), Val{true})
+end
+
+@testset "throw for non positive definite matrix" begin
+    for T in (Float32, Float64, Complex64, Complex128)
+        A = T[1 2; 2 1]; B = T[1, 1]
+        C = cholfact(A)
+        @test_throws PosDefException C\B
+        @test_throws PosDefException det(C)
+        @test_throws PosDefException logdet(C)
+    end
+end