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.

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 tdma with built-in Julia implementation #23

Merged
merged 1 commit into from Jul 20, 2021
Merged

Replace tdma with built-in Julia implementation #23

merged 1 commit into from Jul 20, 2021

Conversation

charleskawczynski
Copy link
Member

🎉

@charleskawczynski
Copy link
Member Author

I locally ran

function tridiag_solve(nz, b_rhs, a, b, c)
    # Note that `1:end` is zero-based indexing.
    A = LinearAlgebra.Tridiagonal(a[1:end], parent(b), c[0:(end - 1)])
    x_new = off_arr(A \ parent(b_rhs))
    x = deepcopy(b_rhs)
    tridiag_solve_old(nz, x, a, b, c)
    if !(x  x_new)
        error("Not matching")
        @show a
        @show b
        @show c
        @show b_rhs
        @show x
        @show x_new
    else
    end
    return x
end


function tridiag_solve_old(nz, x, a, b, c)
    scratch = pyzeros(length(x))
    scratch[0] = c[0] / b[0]
    x[0] = x[0] / b[0]
    @inbounds for i in xrange(1, nz)
        m = 1.0 / (b[i] - a[i] * scratch[i - 1])
        scratch[i] = c[i] * m
        x[i] = (x[i] - a[i] * x[i - 1]) * m
    end
    # TODO: verify translation
    @inbounds for i in revxrange(nz - 2, 0, -1)
        x[i] = x[i] - scratch[i] * x[i + 1]
    end
    return
end

and found that the error is never hit, so the mse table changes are required due to machine precision errors.

Copy link
Member

@haakon-e haakon-e left a comment

Choose a reason for hiding this comment

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

Looks good.

@charleskawczynski
Copy link
Member Author

bors r+

@bors bors bot merged commit f6d3fad into main Jul 20, 2021
@bors bors bot deleted the ck/tdma branch July 20, 2021 17:24
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

Successfully merging this pull request may close these issues.

None yet

2 participants