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"   begin 
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