This code triggers a weird bug on OSX https://github.com/JuliaLang/julia/issues/12899

12899.jl

type BrokenArrowBidiagonal{T} <: AbstractMatrix{T}
    dv::Vector{T}
    av::Vector{T}
    ev::Vector{T}
end
Base.size(B::BrokenArrowBidiagonal) = (n=length(B.dv); (n, n))
function Base.size(B::BrokenArrowBidiagonal, n::Int)
    if n==1 || n==2
        return length(B.dv)
    else
        throw(ArgumentError("invalid dimension $n"))
    end
end
function Base.full{T}(B::BrokenArrowBidiagonal{T})
    n = size(B, 1)
    k = length(B.av)
    M = zeros(T, n, n)
    for i=1:n
        M[i, i] = B.dv[i]
    end
    for i=1:k
        M[i,k+1] = B.av[i]
    end
    for i=k+1:n-1
        M[i, i+1] = B.ev[i-k]
    end
    M
end

Base.svdfact(B::BrokenArrowBidiagonal) = svdfact(full(B))

type PartialFactorization{T, Tr} <: Factorization{T}
    P :: T
    Q :: T
    B :: AbstractMatrix{Tr}
    β :: Tr
end


function thickrestartbidiag(A, q::AbstractVector, l::Int=6, k::Int=2l,
    j::Int=l;
    maxiter::Int=minimum(size(A)), tol::Real=√eps(), reltol::Real=√eps())

    @assert k>l
    L = build(A, q, k)

    local F
    for i=1:maxiter
        info("Iteration $i")
        #@assert size(L.B) == (k, k)
        F = svdfact(L.B)
        L = truncate!(A, L, F, j)
        extend!(A, L, k)
        isconverged(L, F, l, tol, reltol) && break
    end
    F[:S][1:l], L
end

function isconverged(L::PartialFactorization,
        F::Base.LinAlg.SVD, k::Int, tol::Real, reltol::Real)

    @assert tol ≥ 0

    σ = F[:S][1:k]
    Δσ= L.β*abs(F[:U][end, 1:k])

    δσ = copy(Δσ)

    let
        d = Inf
        for i in eachindex(σ), j=1:i-1
            d = min(d, abs(σ[i] - σ[j]))
        end

        #XXX Commenting this loop out makes the error go away
        for i in eachindex(Δσ)
            println("Ritz value ", i, ": ", σ[i])
        end
    end

    all(δσ[1:k] .< max(tol, reltol*σ[1]))
end

#Hernandez2008
function build{T}(A, q::AbstractVector{T}, k::Int)
    m, n = size(A)
    Tr = typeof(real(one(T)))

    P = Array(T, m, k)
    Q = Array(T, n, k+1)
    αs= Array(Tr, k)
    βs= Array(Tr, k-1)
    Q[:, 1] = q
    local β
    p = Array(T, m)
    for j=1:k
        #p = A*q
        A_mul_B!(p, A, q)
        if j>1
            βs[j-1] = β
            Base.LinAlg.axpy!(-β, sub(P, :, j-1), p) #p -= β*P[:,j-1]
        end
        α = norm(p)
        p[:] = p/α

        αs[j] = α
        P[:, j] = p

        Ac_mul_B!(q, A, p) #q = A'p
        Base.LinAlg.axpy!(-1.0, sub(Q, :, 1:j)*(sub(Q, :, 1:j)'q), q) #orthogonalize

        β = norm(q)
        q[:] = q/β
        Q[:,j+1] = q
    end
    PartialFactorization(P, Q, Bidiagonal(αs, βs, true), β)
end

function truncate!(A, L::PartialFactorization,
        F::Base.LinAlg.SVD, k::Int)

    m = size(L.B, 1)
    @assert size(L.P,2)==m==size(L.Q,2)-1

    F0 = svdfact(L.B)
    ρ = L.β*F0[:U][end,:] #Residuals of singular values

    B = [full(L.B) ρ']
    F = svdfact!(B)

    #Take k smallest triplets
    U = F[:U][:,end:-1:end-k+1]
    Σ = F[:S][end:-1:end-k+1]
    V = F[:V][:,end:-1:end-k+1]

    U = F0[:U]*U
    M = eye(m+1)
    M[1:m, 1:m] = F0[:V]
    M = M * V

    r = zeros(m)
    r[end] = 1
    r = - L.β*(L.B\r)

    M = M[1:m, :] + r*M[m+1,:]

    M2 = zeros(m+1, k+1)
    M2[1:m, 1:k] = M
    M2[1:m, k+1] = -r
    M2[m+1, k+1] = 1
    Q, R = qr(M2)

    #XXX Commenting out this block also makes the error go away
    for i=1:size(R,1)
        if R[i,i] < 0
            for j=i:size(R,1)
                R[i,j] = -R[i,j]
            end
        end
    end

    Q = L.Q*Q
    P = L.P*U

    R = (R[1:k+1,1:k] + R[:,k+1]*Q[end,1:k])
    B = (Diagonal(Σ)*triu(R'))
    f = A*sub(L.Q, :, m+1)
    f -= P*ρ[1:k]
    α = norm(f)
    f[:] = f/α
    P = [P f]

    B = UpperTriangular([B; zeros(1,k) α])
    cc=round(B, 3)
    xdump(cc)   # >Commenting either of these out makes the error go away
    display(cc) # >

    g = A'f - α*Q[:, end]
    β = norm(g)
    PartialFactorization(P, Q, B, β)
end

#Hernandez2008
function extend!(A, L::PartialFactorization, k::Int)
    l = size(L.P, 2)-1
    p = L.P[:, end]

    local β
    m, n = size(A)
    q = zeros(n)
    @assert l+1 == size(L.B, 1) == size(L.B, 2)

    if !isa(L.B, Bidiagonal) && !isa(L.B, BrokenArrowBidiagonal) #Cannot be appended
        L.B = [L.B zeros(l+1, k-l-1); zeros(k-l-1, k)]
        @assert size(L.B) == (k, k)
    end

    for j=l+1:k
        Ac_mul_B!(q, A, p) #q = A'p
        q -= L.Q*(L.Q'q)   #orthogonalize
        β = norm(q)
        q[:] = q/β

        L.Q = [L.Q q]
        j==k && break

        A_mul_B!(p, A, q)
        Base.LinAlg.axpy!(-β, sub(L.P, :, j), p)

        α = norm(p)
        p[:] = p/α
        if isa(L.B, Bidiagonal) || isa(L.B, BrokenArrowBidiagonal)
            push!(L.B.dv, α)
            push!(L.B.ev, β)
        else
            L.B[j+1, j+1] = α
            L.B[j  , j+1] = β
        end
        L.P = [L.P p]
    end
    L.β = β
    L
end

let
    srand(1)
    m = 300
    n = 200
    k = 5
    l = 10

    A = randn(m,n)
    q = randn(n)|>x->x/norm(x)
    σ, L = thickrestartbidiag(A, q, k, l, tol=1e-5, maxiter=200)
end