Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Replace magic constant test for linear solve on triangular matrices #5726

Merged
merged 3 commits into from
Feb 8, 2014
Merged

Replace magic constant test for linear solve on triangular matrices #5726

Show file tree
Hide file tree
Changes from 1 commit
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
c82701f
Provides Skeel condition numbers (condskeel)
jiahao Feb 8, 2014
7977571
Replace tests for linear solver on triangular systems
jiahao Feb 8, 2014
76df3e0
Update NEWS.md with condskeel
jiahao Feb 8, 2014
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
Prev Previous commit
Next Next commit
Replace tests for linear solver on triangular systems 
- Add automatic promotion of integers to floats for condskeel
- Ref. #5605, #5705
  • Loading branch information
@jiahao
jiahao committed Feb 8, 2014
commit 797757103ee2b6d235763de76371335a10558088
2 changes: 2 additions & 0 deletions base/linalg/generic.jl
View file
Open in desktop
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)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@jiahao, why were the Integer versions needed, since inv(A) already converts to float?

Note also that abs(inv(A))*abs(A)*abs(x) should be abs(inv(A))*(abs(A)*abs(x)) in order to perform matvec and not matmul operations.

See also the discussion in #17302.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

cc @Sacha0


function issym(A::AbstractMatrix)
m, n = size(A)
Expand Down
32 changes: 29 additions & 3 deletions test/linalg.jl
View file
Open in desktop
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