Skip to content

Commit

Permalink
Replace tests for linear solver on triangular systems
Browse files Browse the repository at this point in the history
- Add automatic promotion of integers to floats for condskeel
- Ref. #5605, #5705
  • Loading branch information
jiahao committed Feb 8, 2014
1 parent c82701f commit 7977571
Show file tree
Hide file tree
Showing 2 changed files with 31 additions and 3 deletions.
2 changes: 2 additions & 0 deletions base/linalg/generic.jl
Original file line number Diff line number Diff line change
Expand Up @@ -176,7 +176,9 @@ cond(x::Number, p) = cond(x)

#Skeel condition numbers
condskeel(A::AbstractMatrix, p::Real=Inf) = norm(abs(inv(A))*abs(A), p)
condskeel{T<:Integer}(A::AbstractMatrix{T}, p::Real=Inf) = norm(abs(inv(float(A)))*abs(A), p)
condskeel(A::AbstractMatrix, x::AbstractVector, p::Real=Inf) = norm(abs(inv(A))*abs(A)*abs(x), p)
condskeel{T<:Integer}(A::AbstractMatrix{T}, x::AbstractVector, p::Real=Inf) = norm(abs(inv(float(A)))*abs(A)*abs(x), p)

function issym(A::AbstractMatrix)
m, n = size(A)
Expand Down
32 changes: 29 additions & 3 deletions test/linalg.jl
Original file line number Diff line number Diff line change
Expand Up @@ -191,13 +191,39 @@ debug && println("Solve square general system of equations")
x = a \ b
@test_approx_eq_eps a*x b 80ε

debug && println("Solve upper trianguler system")
debug && println("Solve upper triangular system")
x = triu(a) \ b
@test_approx_eq_eps triu(a)*x b 20000ε

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

debug && println("Solve lower triangular system")
x = tril(a)\b
@test_approx_eq_eps tril(a)*x b 10000ε

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

debug && println("Test null")
if eltya != BigFloat && eltyb != BigFloat # Revisit when implemented in julia
Expand Down

0 comments on commit 7977571

Please sign in to comment.