tug: julia: Replace sparse matrix with three vectors

src/benchmarks/tug/julia/src/tug/core/btcs.jl

 using FLoops
 using LinearAlgebra
-using SparseArrays
+
+mutable struct Diagonals{T<:Real}
+    left::Vector{T}
+    center::Vector{T}
+    right::Vector{T}
+
+    function Diagonals{T}(size::Int) where {T<:Real}
+        new(Vector{T}(undef, size), Vector{T}(undef, size), Vector{T}(undef, size))
+    end
+end
 
 # calculates coefficient for boundary in constant case
 function calcBoundaryCoeffConstant(
@@ -34,9 +43,9 @@ function createCoeffMatrix(
     numCols::Int,
     rowIndex::Int,
     sx::T,
-)::AbstractMatrix{T} where {T}
+)::Diagonals{T} where {T}
     type::BC_TYPE = BC_TYPE_CLOSED
-    cm::SparseMatrixCSC{T} = spzeros(T, numCols, numCols)
+    diag::Diagonals{T} = Diagonals{T}(numCols)
 
     # left column
     type = getBoundaryElementType(bcLeft[rowIndex])
@@ -45,14 +54,14 @@ function createCoeffMatrix(
         centerCoeffTop, rightCoeffTop =
             calcBoundaryCoeffConstant(alpha[1, rowIndex], alpha[2, rowIndex], sx)
 
-        cm[1, 1] = centerCoeffTop
-        cm[1, 2] = rightCoeffTop
+        diag.center[1] = centerCoeffTop
+        diag.right[1] = rightCoeffTop
     elseif type == BC_TYPE_CLOSED
         centerCoeffTop, rightCoeffTop =
             calcBoundaryCoeffClosed(alpha[1, rowIndex], alpha[2, rowIndex], sx)
 
-        cm[1, 1] = centerCoeffTop
-        cm[1, 2] = rightCoeffTop
+        diag.center[1] = centerCoeffTop
+        diag.right[1] = rightCoeffTop
     else
         error("Undefined Boundary Condition Type somewhere on Left or Top!")
     end
@@ -61,14 +70,14 @@ function createCoeffMatrix(
     n::Int = numCols - 1
 
     for i = 2:n
-        cm[i, i-1] = -sx * calcAlphaIntercell(alpha[i-1, rowIndex], alpha[i, rowIndex])
-        cm[i, i] =
+        diag.left[i] = -sx * calcAlphaIntercell(alpha[i-1, rowIndex], alpha[i, rowIndex])
+        diag.center[i] =
             1 +
             sx * (
                 calcAlphaIntercell(alpha[i, rowIndex], alpha[i+1, rowIndex]) +
                 calcAlphaIntercell(alpha[i-1, rowIndex], alpha[i, rowIndex])
             )
-        cm[i, i+1] = -sx * calcAlphaIntercell(alpha[i, rowIndex], alpha[i+1, rowIndex])
+        diag.right[i] = -sx * calcAlphaIntercell(alpha[i, rowIndex], alpha[i+1, rowIndex])
     end
 
     # right column
@@ -78,19 +87,19 @@ function createCoeffMatrix(
         centerCoeffBottom, leftCoeffBottom =
             calcBoundaryCoeffConstant(alpha[numCols, rowIndex], alpha[numCols-1, rowIndex], sx)
 
-        cm[numCols, numCols-1] = leftCoeffBottom
-        cm[numCols, numCols] = centerCoeffBottom
+        diag.left[numCols] = leftCoeffBottom
+        diag.center[numCols] = centerCoeffBottom
     elseif type == BC_TYPE_CLOSED
         centerCoeffBottom, leftCoeffBottom =
             calcBoundaryCoeffClosed(alpha[numCols, rowIndex], alpha[numCols-1, rowIndex], sx)
 
-        cm[numCols, numCols-1] = leftCoeffBottom
-        cm[numCols, numCols] = centerCoeffBottom
+        diag.left[numCols] = leftCoeffBottom
+        diag.center[numCols] = centerCoeffBottom
     else
         error("Undefined Boundary Condition Type somewhere on Right or Bottom!")
     end
 
-    return cm
+    return diag
 end
 
 # calculates explicit concentration at boundary in closed case
@@ -217,44 +226,24 @@ end
 
 # solver for linear equation system; A corresponds to coefficient matrix, b to the solution vector
 # implementation of Thomas Algorithm
-function ThomasAlgorithm(
-    A::AbstractMatrix{T},
-    x_vec::AbstractVector{T},
-)::AbstractVector{T} where {T}
+function ThomasAlgorithm(A::Diagonals{T}, x_vec::AbstractVector{T})::AbstractVector{T} where {T}
     n::Int = size(x_vec, 1)
 
-    a_diag = Vector{T}(undef, n)
-    b_diag = Vector{T}(undef, n)
-    c_diag = Vector{T}(undef, n)
-
-    # Fill diagonals vectors
-    b_diag[1] = A[1, 1]
-    c_diag[1] = A[1, 2]
-
-    for i = 2:n-1
-        a_diag[i] = A[i, i-1]
-        b_diag[i] = A[i, i]
-        c_diag[i] = A[i, i+1]
-    end
-
-    a_diag[n] = A[n, n-1]
-    b_diag[n] = A[n, n]
-
     # start solving - c_diag and x_vec are overwritten
     n -= 1
 
-    c_diag[1] = c_diag[1] / b_diag[1]
-    x_vec[1] = x_vec[1] / b_diag[1]
+    A.right[1] = A.right[1] / A.center[1]
+    x_vec[1] = x_vec[1] / A.center[1]
 
     for i = 2:n
-        c_diag[i] /= b_diag[i] - a_diag[i] * c_diag[i-1]
-        x_vec[i] = (x_vec[i] - a_diag[i] * x_vec[i-1]) / (b_diag[i] - a_diag[i] * c_diag[i-1])
+        A.right[i] /= A.center[i] - A.left[i] * A.right[i-1]
+        x_vec[i] = (x_vec[i] - A.left[i] * x_vec[i-1]) / (A.center[i] - A.left[i] * A.right[i-1])
     end
 
-    x_vec[n+1] = (x_vec[n+1] - a_diag[n+1] * x_vec[n]) / (b_diag[n+1] - a_diag[n+1] * c_diag[n])
+    x_vec[n+1] = (x_vec[n+1] - A.left[n+1] * x_vec[n]) / (A.center[n+1] - A.left[n+1] * A.right[n])
 
     for i = 0:n-1
-        x_vec[n-i] -= c_diag[n-i] * x_vec[n-i+1]
+        x_vec[n-i] -= A.right[n-i] * x_vec[n-i+1]
     end
 
     return x_vec