X = rand(5,5)
A = X + X'
display(A)
for n = 0:10
    Q,R = qr(A)
    A = R*Q
    display(A)
end

eigvals(X + X')

function matadd1!(C, A, B)
    m,n = size(A)
    (m,n) == size(B) || error("matrix sizes must match")
    for i = 1:m
        for j = 1:n
            @inbounds C[i,j] = A[i,j] + B[i,j]
        end
    end
    return C
end
matadd1(A,B) = matadd1!(similar(A, promote_type(eltype(A),eltype(B))), A,B)

# same as matadd1 except the loop order is reversed
function matadd2!(C, A, B)
    m,n = size(A)
    (m,n) == size(B) || error("matrix sizes must match")
    for j = 1:n
        for i = 1:m
            @inbounds C[i,j] = A[i,j] + B[i,j]
        end
    end
    return C
end
matadd2(A,B) = matadd2!(similar(A, promote_type(eltype(A),eltype(B))), A,B)

# correctness check:
A = rand(5,7)
B = rand(5,7)
norm(matadd1(A,B) - (A+B)), norm(matadd2(A,B) - (A+B))

?A+B

N = round(Int, linspace(10, 1000, 200))
t0 = Float64[]
t1 = Float64[]
t2 = Float64[]
for n in N
    A = rand(n,n)
    B = rand(n,n)
    C = similar(A)
    push!(t0, @elapsed A+B) # note: does not preallocate output
    push!(t1, @elapsed matadd1!(C,A,B))
    push!(t2, @elapsed matadd2!(C,A,B))
end

using PyPlot

plot(N, N.^2 ./ t0 * 1e-9, "k.-")
plot(N, N.^2 ./ t1 * 1e-9, "bo-")
plot(N, N.^2 ./ t2 * 1e-9, "rs-")
xlabel(L"n")
ylabel(L"gflops ($n^2/t$)")
legend(["A+B","matmul1!","matmul2!"])