Skip to content

Commit

Permalink
Remove old debug functionality in triangular.jl (#1124)
Browse files Browse the repository at this point in the history 
Co-authored-by: Daniel Karrasch <daniel.karrasch@posteo.de>
  • Loading branch information
@andreasnoack @dkarrasch
andreasnoack and dkarrasch authored Dec 22, 2024
1 parent 8d9b14f commit 959d985
Showing 1 changed file with 53 additions and 64 deletions.
117 changes: 53 additions & 64 deletions test/triangular.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,6 @@

module TestTriangular

debug = false
using Test, LinearAlgebra, Random
using LinearAlgebra: BlasFloat, errorbounds, full!, transpose!,
UnitUpperTriangular, UnitLowerTriangular,
Expand All @@ -16,14 +15,13 @@ using .Main.SizedArrays
isdefined(Main, :FillArrays) || @eval Main include(joinpath($(BASE_TEST_PATH), "testhelpers", "FillArrays.jl"))
using .Main.FillArrays

debug && println("Triangular matrices")

n = 9
Random.seed!(123)

debug && println("Test basic type functionality")
@test_throws DimensionMismatch LowerTriangular(randn(5, 4))
@test LowerTriangular(randn(3, 3)) |> t -> [size(t, i) for i = 1:3] == [size(Matrix(t), i) for i = 1:3]
@testset "Test basic type functionality" begin
@test_throws DimensionMismatch LowerTriangular(randn(5, 4))
@test LowerTriangular(randn(3, 3)) |> t -> [size(t, i) for i = 1:3] == [size(Matrix(t), i) for i = 1:3]
end

struct MyTriangular{T, A<:LinearAlgebra.AbstractTriangular{T}} <: LinearAlgebra.AbstractTriangular{T}
data :: A
Expand Down Expand Up @@ -51,8 +49,6 @@ Base.getindex(A::MyTriangular, i::Int, j::Int) = A.data[i,j]
t1t = t1{elty1}(t1(rand(Int8, n, n)))
@test typeof(t1t) == t1{elty1, Matrix{elty1}}

debug && println("elty1: $elty1, A1: $t1")

# Convert
@test convert(AbstractMatrix{elty1}, A1) == A1
@test convert(Matrix, A1) == A1
Expand Down Expand Up @@ -356,8 +352,6 @@ Base.getindex(A::MyTriangular, i::Int, j::Int) = A.data[i,j]
(LowerTriangular, :L),
(UnitLowerTriangular, :L))

debug && println("elty1: $elty1, A1: $t1, elty2: $elty2, A2: $t2")

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')))
M2 = Matrix(A2)
# Convert
Expand Down Expand Up @@ -451,8 +445,6 @@ Base.getindex(A::MyTriangular, i::Int, j::Int) = A.data[i,j]

B = convert(Matrix{eltyB}, (elty1 <: Complex ? real(A1) : A1)*fill(1., n, n))

debug && println("elty1: $elty1, A1: $t1, B: $eltyB")

Tri = Tridiagonal(rand(eltyB,n-1),rand(eltyB,n),rand(eltyB,n-1))
C = Matrix{promote_type(elty1,eltyB)}(undef, n, n)
mul!(C, Tri, A1)
Expand Down Expand Up @@ -639,83 +631,80 @@ Aimg = randn(n, n)/2
A2real = randn(n, n)/2
A2img = randn(n, n)/2

for eltya in (Float32, Float64, ComplexF32, ComplexF64, BigFloat, Int)
@testset 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))))

for eltyb in (Float32, Float64, ComplexF32, ComplexF64)
εb = eps(abs(float(one(eltyb))))
ε = max(εa,εb)

debug && println("\ntype of A: ", eltya, " type of b: ", eltyb, "\n")
@testset "Solve upper triangular system" begin
Atri = UpperTriangular(lu(A).U) |> t -> eltya <: Complex && eltyb <: Real ? real(t) : t # Here the triangular matrix can't be too badly conditioned
b = convert(Matrix{eltyb}, Matrix(Atri)*fill(1., n, 2))
x = Atri \ b

debug && println("Solve upper triangular system")
Atri = UpperTriangular(lu(A).U) |> t -> eltya <: Complex && eltyb <: Real ? real(t) : t # Here the triangular matrix can't be too badly conditioned
b = convert(Matrix{eltyb}, Matrix(Atri)*fill(1., n, 2))
x = Matrix(Atri) \ b
# Test error estimates
if eltya != BigFloat && eltyb != BigFloat
for i = 1:2
@test norm(x[:,1] .- 1) <= errorbounds(UpperTriangular(A), x, b)[1][i]
end
end

debug && println("Test error estimates")
if eltya != BigFloat && eltyb != BigFloat
for i = 1:2
@test norm(x[:,1] .- 1) <= errorbounds(UpperTriangular(A), x, b)[1][i]
# Test forward error [JIN 5705] if this is not a BigFloat
γ = n*ε/(1 - n*ε)
if eltya != BigFloat
bigA = big.(Atri)
x̂ = fill(1., n, 2)
for i = 1:size(b, 2)
@test norm(x̂[:,i] - x[:,i], Inf)/norm(x̂[:,i], Inf) <= condskeel(bigA, x̂[:,i])*γ/(1 - condskeel(bigA)*γ)
end
end
end
debug && println("Test forward error [JIN 5705] if this is not a BigFloat")

x = Atri \ b
γ = n*ε/(1 - n*ε)
if eltya != BigFloat
bigA = big.(Atri)
x̂ = fill(1., n, 2)
# Test backward error [JIN 5705]
for i = 1:size(b, 2)
@test norm(x̂[:,i] - x[:,i], Inf)/norm(x̂[:,i], Inf) <= condskeel(bigA, x̂[:,i])*γ/(1 - condskeel(bigA)*γ)
@test norm(abs.(b[:,i] - Atri*x[:,i]), Inf) <= γ * norm(Atri, Inf) * norm(x[:,i], Inf)
end
end

debug && println("Test backward error [JIN 5705]")
for i = 1:size(b, 2)
@test norm(abs.(b[:,i] - Atri*x[:,i]), Inf) <= γ * norm(Atri, Inf) * norm(x[:,i], Inf)
end
@testset "Solve lower triangular system" begin
Atri = LowerTriangular(lu(A).L) |> t -> eltya <: Complex && eltyb <: Real ? real(t) : t # Here the triangular matrix can't be too badly conditioned
b = convert(Matrix{eltyb}, Matrix(Atri)*fill(1., n, 2))
x = Atri \ b

debug && println("Solve lower triangular system")
Atri = UpperTriangular(lu(A).U) |> t -> eltya <: Complex && eltyb <: Real ? real(t) : t # Here the triangular matrix can't be too badly conditioned
b = convert(Matrix{eltyb}, Matrix(Atri)*fill(1., n, 2))
x = Matrix(Atri)\b
# Test error estimates
if eltya != BigFloat && eltyb != BigFloat
for i = 1:2
@test norm(x[:,1] .- 1) <= errorbounds(LowerTriangular(A), x, b)[1][i]
end
end

debug && println("Test error estimates")
if eltya != BigFloat && eltyb != BigFloat
for i = 1:2
@test norm(x[:,1] .- 1) <= errorbounds(UpperTriangular(A), x, b)[1][i]
# Test forward error [JIN 5705] if this is not a BigFloat
γ = n*ε/(1 - n*ε)
if eltya != BigFloat
bigA = big.(Atri)
x̂ = fill(1., n, 2)
for i = 1:size(b, 2)
@test norm(x̂[:,i] - x[:,i], Inf)/norm(x̂[:,i], Inf) <= condskeel(bigA, x̂[:,i])*γ/(1 - condskeel(bigA)*γ)
end
end
end

debug && println("Test forward error [JIN 5705] if this is not a BigFloat")
b = (b0 = Atri*fill(1, n, 2); convert(Matrix{eltyb}, eltyb == Int ? trunc.(b0) : b0))
x = Atri \ b
γ = n*ε/(1 - n*ε)
if eltya != BigFloat
bigA = big.(Atri)
x̂ = fill(1., n, 2)
# Test backward error [JIN 5705]
for i = 1:size(b, 2)
@test norm(x̂[:,i] - x[:,i], Inf)/norm(x̂[:,i], Inf) <= condskeel(bigA, x̂[:,i])*γ/(1 - condskeel(bigA)*γ)
@test norm(abs.(b[:,i] - Atri*x[:,i]), Inf) <= γ * norm(Atri, Inf) * norm(x[:,i], Inf)
end
end

debug && println("Test backward error [JIN 5705]")
for i = 1:size(b, 2)
@test norm(abs.(b[:,i] - Atri*x[:,i]), Inf) <= γ * norm(Atri, Inf) * norm(x[:,i], Inf)
end
end
end

# Issue 10742 and similar
@test istril(UpperTriangular(diagm(0 => [1,2,3,4])))
@test istriu(LowerTriangular(diagm(0 => [1,2,3,4])))
@test isdiag(UpperTriangular(diagm(0 => [1,2,3,4])))
@test isdiag(LowerTriangular(diagm(0 => [1,2,3,4])))
@test !isdiag(UpperTriangular(rand(4, 4)))
@test !isdiag(LowerTriangular(rand(4, 4)))
@testset "triangularity/diagonality of triangular views (#10742)" begin
@test istril(UpperTriangular(diagm(0 => [1,2,3,4])))
@test istriu(LowerTriangular(diagm(0 => [1,2,3,4])))
@test isdiag(UpperTriangular(diagm(0 => [1,2,3,4])))
@test isdiag(LowerTriangular(diagm(0 => [1,2,3,4])))
@test !isdiag(UpperTriangular(rand(4, 4)))
@test !isdiag(LowerTriangular(rand(4, 4)))
end

# Test throwing in fallbacks for non BlasFloat/BlasComplex in A_rdiv_Bx!
let n = 5
Expand Down

0 comments on commit 959d985

Please sign in to comment.