In [1]:
## Lecture 6 - Simple linear system solving!

"""
`diag_solve(A,b) -> x`

Solve a linear system of equations with a diagonal matrix A.

`x = diag_solve(A,b)` returns a vector x such that Ax=b when
A is a diagonal n-by-n matrix. If A is not diagonal, then this
routine may silently fail to return such an x; use with
extreme care!

Example
-------
    n = 6
    A = diagm(randn(n))
    b = randn(n)
    x = diag_solve(A,b)
    norm(x - A\b) # should be close to zero!
"""
function diag_solve(A,b)
    n = size(A,1)
    @assert n == length(b) "b has the wrong length"
    @assert n == size(A,2) "A is not square"

    x = zeros(n) # allocate the solution
    for i=1:n
        @assert A[i,i] > eps(1.0) "A is singular"
        if abs(A[i,i]) <= eps(1.0)
            error("A is singular")
        end
        x[i] = b[i]/A[i,i]
    end
    return x
end
Out[1]:
diag_solve
In [2]:
##
n = 6
A = diagm(randn(n))
b = randn(n)
x = diag_solve(A,b)
norm(x - A\b) # should be close to zero!
AssertionError: A is singular

Stacktrace:
 [1] diag_solve(::Array{Float64,2}, ::Array{Float64,1}) at ./In[1]:28
 [2] include_string(::String, ::String) at ./loading.jl:515
In [3]:
## is this fast? Timing version 1
n = 50
A = diagm(randn(n))
b = randn(n)
@time diag_solve(A,b)
@time A\b
AssertionError: A is singular

Stacktrace:
 [1] diag_solve(::Array{Float64,2}, ::Array{Float64,1}) at ./In[1]:28
 [2] include_string(::String, ::String) at ./loading.jl:515
In [4]:
## Timing version 2

""" Run a time trial on solving n-by-n systems, t times. """
function many_times(n,t)
    time_us = 0.0
    time_julia = 0.0

    for i=1:t
        A = diagm(randn(n))
        b = randn(n)
        time_us += @elapsed diag_solve(A,b)
        time_julia += @elapsed A\b
    end
    return time_us, time_julia
end
Out[4]:
many_times
In [5]:
##
@show many_times(1000,1000)
AssertionError: A is singular

Stacktrace:
 [1] diag_solve(::Array{Float64,2}, ::Array{Float64,1}) at ./In[1]:28
 [2] macro expansion at ./util.jl:293 [inlined]
 [3] many_times(::Int64, ::Int64) at ./In[4]:11
 [4] include_string(::String, ::String) at ./loading.jl:515
In [6]:
## Does backslash work with sparse matrices.
function many_times_sparse(n,t)
    time_us = 0.0
    time_julia = 0.0

    for i=1:t
        A = diagm(randn(n))
        b = randn(n)
        As = sparse(A) # convert into a sparse matrix
        time_us += @elapsed diag_solve(As,b)
        time_julia += @elapsed As\b
    end
    return time_us, time_julia
end
@show many_times_sparse(100,1000)
AssertionError: A is singular

Stacktrace:
 [1] diag_solve(::SparseMatrixCSC{Float64,Int64}, ::Array{Float64,1}) at ./In[1]:28
 [2] macro expansion at ./util.jl:293 [inlined]
 [3] many_times_sparse(::Int64, ::Int64) at ./In[6]:10
 [4] include_string(::String, ::String) at ./loading.jl:515