Fixing Julia for Thomas algorithm (#495)

Author: leios

contents/thomas_algorithm/code/julia/thomas.jl

@@ -1,23 +1,45 @@
-function(a::Vector{Float64}, b::Vector{Float64}, c::Vector{Float64},
-         d::Vector{Float64}, soln::Vector{Float64})
+function thomas(a::Vector{Float64}, b::Vector{Float64}, c::Vector{Float64},
+                d::Vector{Float64}, n::Int64)
+
+    x = copy(d)
+    c_prime = copy(c)
+
     # Setting initial elements
-    c[0] = c[0] / b[0]
-    d[0] = d[0] / b[0]
-
-    for i = 1:n
-        # Scale factor is for c and d
-        scale = 1 / (b[i] - c[i-1]*a[i])
-        c[i] = c[i] * scale
-        d[i] = (d[i] - a[i] * d[i-1]) * scale
-    end
+    c_prime[1] /= b[1]
+    x[1] /= b[1]
 
-    # Set the last solution for back-substitution
-    soln[n-1] = d[n-1]
+    for i = 2:n
+        # Scale factor is for c_prime and x
+        scale = 1.0 / (b[i] - c_prime[i-1]*a[i])
+        c_prime[i] *= scale
+        x[i] = (x[i] - a[i] * x[i-1]) * scale
+    end
 
     # Back-substitution
-    for i = n-2:0
-        soln[i] = d[i] - c[i] * soln[i+1]
+    for i = n-1:-1:1
+        x[i] -= (c_prime[i] * x[i+1])
     end
 
+    return x
+
+end
+
+function main()
+    a = [0.0, 2.0, 3.0]
+    b = [1.0, 3.0, 6.0]
+    c = [4.0, 5.0, 0.0]
+    d = [7.0, 5.0, 3.0]
+
+    println("The system,")
+    println(a)
+    println(b)
+    println(c)
+    println(d)
+    println("Has the solution:")
+
+    soln = thomas(a, b, c, d, 3)
+
+    println(soln)
 end
 
+main()