22
33module TestTriangular
44
5- debug = false
65using Test, LinearAlgebra, Random
76using LinearAlgebra: BlasFloat, errorbounds, full!, transpose!,
87 UnitUpperTriangular, UnitLowerTriangular,
@@ -16,14 +15,13 @@ using .Main.SizedArrays
1615isdefined (Main, :FillArrays ) || @eval Main include (joinpath ($ (BASE_TEST_PATH), " testhelpers" , " FillArrays.jl" ))
1716using . Main. FillArrays
1817
19- debug && println (" Triangular matrices" )
20-
2118n = 9
2219Random. seed! (123 )
2320
24- debug && println (" Test basic type functionality" )
25- @test_throws DimensionMismatch LowerTriangular (randn (5 , 4 ))
26- @test LowerTriangular (randn (3 , 3 )) |> t -> [size (t, i) for i = 1 : 3 ] == [size (Matrix (t), i) for i = 1 : 3 ]
21+ @testset " Test basic type functionality"
22+ @test_throws DimensionMismatch LowerTriangular (randn (5 , 4 ))
23+ @test LowerTriangular (randn (3 , 3 )) |> t -> [size (t, i) for i = 1 : 3 ] == [size (Matrix (t), i) for i = 1 : 3 ]
24+ end
2725
2826struct MyTriangular{T, A<: LinearAlgebra.AbstractTriangular{T} } <: LinearAlgebra.AbstractTriangular{T}
2927 data :: A
@@ -51,8 +49,6 @@ Base.getindex(A::MyTriangular, i::Int, j::Int) = A.data[i,j]
5149 t1t = t1 {elty1} (t1 (rand (Int8, n, n)))
5250 @test typeof (t1t) == t1{elty1, Matrix{elty1}}
5351
54- debug && println (" elty1: $elty1 , A1: $t1 " )
55-
5652 # Convert
5753 @test convert (AbstractMatrix{elty1}, A1) == A1
5854 @test convert (Matrix, A1) == A1
@@ -351,8 +347,6 @@ Base.getindex(A::MyTriangular, i::Int, j::Int) = A.data[i,j]
351347 (LowerTriangular, :L ),
352348 (UnitLowerTriangular, :L ))
353349
354- debug && println (" elty1: $elty1 , A1: $t1 , elty2: $elty2 , A2: $t2 " )
355-
356350 A2 = t2 (elty2 == Int ? rand (1 : 7 , n, n) : convert (Matrix{elty2}, (elty2 <: Complex ? complex .(randn (n, n), randn (n, n)) : randn (n, n)) |> t -> cholesky (t' t). U |> t -> uplo2 === :U ? t : copy (t' )))
357351 M2 = Matrix (A2)
358352 # Convert
@@ -440,8 +434,6 @@ Base.getindex(A::MyTriangular, i::Int, j::Int) = A.data[i,j]
440434 for eltyB in (Float32, Float64, BigFloat, ComplexF32, ComplexF64, Complex{BigFloat})
441435 B = convert (Matrix{eltyB}, (elty1 <: Complex ? real (A1) : A1)* fill (1. , n, n))
442436
443- debug && println (" elty1: $elty1 , A1: $t1 , B: $eltyB " )
444-
445437 Tri = Tridiagonal (rand (eltyB,n- 1 ),rand (eltyB,n),rand (eltyB,n- 1 ))
446438 C = Matrix {promote_type(elty1,eltyB)} (undef, n, n)
447439 mul! (C, Tri, A1)
@@ -628,73 +620,69 @@ Aimg = randn(n, n)/2
628620A2real = randn (n, n)/ 2
629621A2img = randn (n, n)/ 2
630622
631- for eltya in (Float32, Float64, ComplexF32, ComplexF64, BigFloat, Int)
623+ @testset for eltya in (Float32, Float64, ComplexF32, ComplexF64, BigFloat, Int)
632624 A = eltya == Int ? rand (1 : 7 , n, n) : convert (Matrix{eltya}, eltya <: Complex ? complex .(Areal, Aimg) : Areal)
633- # a2 = eltya == Int ? rand(1:7, n, n) : convert(Matrix{eltya}, eltya <: Complex ? complex.(a2real, a2img) : a2real)
634625 εa = eps (abs (float (one (eltya))))
635626
636627 for eltyb in (Float32, Float64, ComplexF32, ComplexF64)
637628 εb = eps (abs (float (one (eltyb))))
638629 ε = max (εa,εb)
639630
640- debug && println (" \n type of A: " , eltya, " type of b: " , eltyb, " \n " )
631+ @testset " Solve upper triangular system" begin
632+ Atri = UpperTriangular (lu (A). U) |> t -> eltya <: Complex && eltyb <: Real ? real (t) : t # Here the triangular matrix can't be too badly conditioned
633+ b = convert (Matrix{eltyb}, Matrix (Atri)* fill (1. , n, 2 ))
634+ x = Atri \ b
641635
642- debug && println (" Solve upper triangular system" )
643- Atri = UpperTriangular (lu (A). U) |> t -> eltya <: Complex && eltyb <: Real ? real (t) : t # Here the triangular matrix can't be too badly conditioned
644- b = convert (Matrix{eltyb}, Matrix (Atri)* fill (1. , n, 2 ))
645- x = Matrix (Atri) \ b
636+ # Test error estimates
637+ if eltya != BigFloat && eltyb != BigFloat
638+ for i = 1 : 2
639+ @test norm (x[:,1 ] .- 1 ) <= errorbounds (UpperTriangular (A), x, b)[1 ][i]
640+ end
641+ end
646642
647- debug && println (" Test error estimates" )
648- if eltya != BigFloat && eltyb != BigFloat
649- for i = 1 : 2
650- @test norm (x[:,1 ] .- 1 ) <= errorbounds (UpperTriangular (A), x, b)[1 ][i]
643+ # Test forward error [JIN 5705] if this is not a BigFloat
644+ γ = n* ε/ (1 - n* ε)
645+ if eltya != BigFloat
646+ bigA = big .(Atri)
647+ x̂ = fill (1. , n, 2 )
648+ for i = 1 : size (b, 2 )
649+ @test norm (x̂[:,i] - x[:,i], Inf )/ norm (x̂[:,i], Inf ) <= condskeel (bigA, x̂[:,i])* γ/ (1 - condskeel (bigA)* γ)
650+ end
651651 end
652- end
653- debug && println (" Test forward error [JIN 5705] if this is not a BigFloat" )
654652
655- x = Atri \ b
656- γ = n* ε/ (1 - n* ε)
657- if eltya != BigFloat
658- bigA = big .(Atri)
659- x̂ = fill (1. , n, 2 )
653+ # Test backward error [JIN 5705]
660654 for i = 1 : size (b, 2 )
661- @test norm (x̂ [:,i] - x[:,i], Inf )/ norm (x̂[:,i], Inf ) <= condskeel (bigA, x̂ [:,i]) * γ / ( 1 - condskeel (bigA) * γ )
655+ @test norm (abs .(b [:,i] - Atri * x[:,i]) , Inf ) <= γ * norm (Atri, Inf ) * norm (x [:,i], Inf )
662656 end
663657 end
664658
665- debug && println ( " Test backward error [JIN 5705] " )
666- for i = 1 : size (b, 2 )
667- @test norm ( abs .(b[:,i] - Atri * x[:,i]), Inf ) <= γ * norm (Atri, Inf ) * norm (x[:,i], Inf )
668- end
659+ @testset " Solve lower triangular system " begin
660+ Atri = UpperTriangular ( lu (A) . U) |> t -> eltya <: Complex && eltyb <: Real ? real (t) : t # Here the triangular matrix can't be too badly conditioned
661+ b = convert (Matrix{eltyb}, Matrix (Atri) * fill ( 1. , n, 2 ) )
662+ x = Atri \ b
669663
670- debug && println (" Solve lower triangular system" )
671- Atri = UpperTriangular (lu (A). U) |> t -> eltya <: Complex && eltyb <: Real ? real (t) : t # Here the triangular matrix can't be too badly conditioned
672- b = convert (Matrix{eltyb}, Matrix (Atri)* fill (1. , n, 2 ))
673- x = Matrix (Atri)\ b
664+ # Test error estimates
665+ if eltya != BigFloat && eltyb != BigFloat
666+ for i = 1 : 2
667+ @test norm (x[:,1 ] .- 1 ) <= errorbounds (UpperTriangular (A), x, b)[1 ][i]
668+ end
669+ end
674670
675- debug && println (" Test error estimates" )
676- if eltya != BigFloat && eltyb != BigFloat
677- for i = 1 : 2
678- @test norm (x[:,1 ] .- 1 ) <= errorbounds (UpperTriangular (A), x, b)[1 ][i]
671+ # Test forward error [JIN 5705] if this is not a BigFloat
672+ γ = n* ε/ (1 - n* ε)
673+ if eltya != BigFloat
674+ bigA = big .(Atri)
675+ x̂ = fill (1. , n, 2 )
676+ for i = 1 : size (b, 2 )
677+ @test norm (x̂[:,i] - x[:,i], Inf )/ norm (x̂[:,i], Inf ) <= condskeel (bigA, x̂[:,i])* γ/ (1 - condskeel (bigA)* γ)
678+ end
679679 end
680- end
681680
682- debug && println (" Test forward error [JIN 5705] if this is not a BigFloat" )
683- b = (b0 = Atri* fill (1 , n, 2 ); convert (Matrix{eltyb}, eltyb == Int ? trunc .(b0) : b0))
684- x = Atri \ b
685- γ = n* ε/ (1 - n* ε)
686- if eltya != BigFloat
687- bigA = big .(Atri)
688- x̂ = fill (1. , n, 2 )
681+ # Test backward error [JIN 5705]
689682 for i = 1 : size (b, 2 )
690- @test norm (x̂ [:,i] - x[:,i], Inf )/ norm (x̂[:,i], Inf ) <= condskeel (bigA, x̂ [:,i]) * γ / ( 1 - condskeel (bigA) * γ )
683+ @test norm (abs .(b [:,i] - Atri * x[:,i]) , Inf ) <= γ * norm (Atri, Inf ) * norm (x [:,i], Inf )
691684 end
692685 end
693-
694- debug && println (" Test backward error [JIN 5705]" )
695- for i = 1 : size (b, 2 )
696- @test norm (abs .(b[:,i] - Atri* x[:,i]), Inf ) <= γ * norm (Atri, Inf ) * norm (x[:,i], Inf )
697- end
698686 end
699687end
700688
0 commit comments