tridiagonal.h

Go to the documentation of this file.

#ifndef FUSED_ALGORITHMS_TRIDIAGONAL_H
#define FUSED_ALGORITHMS_TRIDIAGONAL_H
 
#include "config.h"
#include "utilities/expected.h"
#include "matrix_traits.h"
#include "fused/fused_vector.h"
#include "math/math_utils.h"
 
namespace matrix_algorithms
{
 
    using matrix_traits::MatrixStatus;
 
    template <typename T, my_size_t N>
    Expected<FusedVector<T, N>, MatrixStatus> thomas_solve(
        const FusedVector<T, N> &a,
        const FusedVector<T, N> &d,
        const FusedVector<T, N> &c,
        const FusedVector<T, N> &b)
    {
        static_assert(is_floating_point_v<T>,
                      "thomas_solve requires a floating-point scalar type");
 
        // Work on copies (modified during forward sweep)
        FusedVector<T, N> dp = d;
        FusedVector<T, N> bp = b;
 
        // Forward sweep
        for (my_size_t i = 1; i < N; ++i)
        {
            T diag = dp(i - 1);
 
            if (math::abs(diag) <= T(PRECISION_TOLERANCE))
            {
                return Unexpected{MatrixStatus::Singular};
            }
 
            T w = a(i) / diag;
            dp(i) = dp(i) - w * c(i - 1);
            bp(i) = bp(i) - w * bp(i - 1);
        }
 
        // Check last pivot
        if (math::abs(dp(N - 1)) <= T(PRECISION_TOLERANCE))
        {
            return Unexpected{MatrixStatus::Singular};
        }
 
        // Back substitution
        FusedVector<T, N> x(T(0));
 
        x(N - 1) = bp(N - 1) / dp(N - 1);
 
        for (my_size_t i = N - 1; i-- > 0;)
        {
            x(i) = (bp(i) - c(i) * x(i + 1)) / dp(i);
        }
 
        return move(x);
    }
 
} // namespace matrix_algorithms
 
#endif // FUSED_ALGORITHMS_TRIDIAGONAL_H