# Custom linear system solver

**URL:** <https://discourse.julialang.org/t/custom-linear-system-solver/22094>\
**Category:** New to Julia\
**Tags:** question\
**Created:** [March 20, 2019, 8:46am UTC](https://discourse.julialang.org/t/custom-linear-system-solver/22094 "2019-03-20T08:46:16Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![pupuyuyu](https://avatars.discourse-cdn.com/v4/letter/p/5f9b8f/32.png) [@pupuyuyu](https://discourse.julialang.org/u/pupuyuyu)\
**Post date:** [March 20, 2019, 8:46am UTC](https://discourse.julialang.org/t/custom-linear-system-solver/22094/1 "2019-03-20T08:46:16Z")

</div>

Hi there,

I am new to julia. I am trying to write a function to solve tridiagonal matrix system by [Thomas Algorithm](https://en.wikipedia.org/wiki/Tridiagonal_matrix_algorithm#Method)

Here is my code

```julia
function TDMA(a,b,c,D)
# Thomas Algorithm
# Vectors a,b,c have the same length
# Should have a[1]=0 and c[n]=0
# A*x = D
# [b[1] c[1] ] [x[1] ] [D[1] ]
# [a[2] b[2] c[2] ] [x[2] ] [D[2] ]
# [... ... ...] [...] = [...]
# [a[n-1] b[n-1] c[n-1]] [x[n-1]] [D[n-1]]
# [a[n] b[n] ] [x[n] ] [D[n] ]  
x = zeros(length(D),1)
d = D
n = length(D)
w = zeros(n,1)
    for i = 2:n
        w[i] = a[i]/b[i-1]
        b[i] = b[i] - w[i]*c[i-1]
        d[i] = d[i] - w[i]*d[i-1]
    end
    x[n] = d[n]/b[n]
    for k = 1:n-1
        x[n-k] = (d[n-k]-c[n-k]*x[n-k+1])/b[n-k]
    end
    x
end

```

I am testing it by this example

```julia
A = 
[1 2 0;
 4 5 6;
 0 1 2]
Y = [1; 3; 4]

```

with

```julia
a = [0 4 1]
b = [1 5 2]
c = [2 6 0]

```

to be the arguments.

The build-in function

```julia
A\Y

```

gives the solution

```julia
 -3.333333333333333 
  2.1666666666666665
  0.9166666666666666

```

Howevery, when I try to use my TDMA to solve the equation, it comes error

```julia
TDMA(a,b,c,Y)

InexactError: Int64(4.545454545454545)

Stacktrace:
 [1] Type at ./float.jl:703 [inlined]
 [2] convert at ./number.jl:7 [inlined]
 [3] setindex! at ./array.jl:767 [inlined]
 [4] TDMA(::Array{Int64,2}, ::Array{Int64,2}, ::Array{Int64,2}, ::Array{Int64,1}) at ./In[13]:17
 [5] top-level scope at In[14]:1

```

Can anyone help me with this?

Thanks

---

<div class="post-metadata">

**Author:** ![jbrea](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jbrea/32/3879_2.png) [@jbrea](https://discourse.julialang.org/u/jbrea)\
**Post date:** [March 20, 2019, 9:23am UTC](https://discourse.julialang.org/t/custom-linear-system-solver/22094/2 "2019-03-20T09:23:16Z")

</div>

Hi

The problem is the assignment of a floating point number to a vector of integers in the lines `b[i] = b[i] - w[i]*c[i-1]` and `d[i] = d[i] - w[i]*d[i-1]`.

You see a similar error when you try to do

```julia
julia> b = [1, 2, 3]
3-element Array{Int64,1}:
 1
 2
 3

julia> b[1] = .3
ERROR: InexactError: Int64(0.3)
Stacktrace:
 [1] Type at ./float.jl:703 [inlined]
 [2] convert at ./number.jl:7 [inlined]
 [3] setindex!(::Array{Int64,1}, ::Float64, ::Int64) at ./array.jl:767
 [4] top-level scope at none:0

```

The solution is to either call your function with arrays of floats, e.g. by writing `b = [1., 5., 2.]`, or convert these vectors inside you function with e.g. `b = float.(b)` and `d = float.(D)`.  
Note also that it is common to [append `!` to function names that change the arguments](https://docs.julialang.org/en/v1/manual/style-guide/#Append-!-to-names-of-functions-that-modify-their-arguments-1).

---

<div class="post-metadata">

**Author:** ![pupuyuyu](https://avatars.discourse-cdn.com/v4/letter/p/5f9b8f/32.png) [@pupuyuyu](https://discourse.julialang.org/u/pupuyuyu)\
**Post date:** [March 20, 2019, 5:53pm UTC](https://discourse.julialang.org/t/custom-linear-system-solver/22094/3 "2019-03-20T17:53:24Z")

</div>

It works now. Thank you!

---

<div class="post-metadata">

**Author:** ![tkoolen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tkoolen/32/1603_2.png) [@tkoolen](https://discourse.julialang.org/u/tkoolen)\
**Post date:** [March 20, 2019, 7:57pm UTC](https://discourse.julialang.org/t/custom-linear-system-solver/22094/4 "2019-03-20T19:57:53Z")

</div>

Note that there is a `Tridiagonal` type in LinearAlgebra by the way: [Linear Algebra · The Julia Language](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/#LinearAlgebra.Tridiagonal). `\` is probably overloaded for this type.

---

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [December 9, 2020, 9:02am UTC](https://discourse.julialang.org/t/custom-linear-system-solver/22094/5 "2020-12-09T09:02:46Z")

</div>

Yes but to me, it is slower than thomas. However, I think TDMA is numerically less stable