thomas.cpp

Go to the documentation of this file.

/// @file linalg/factorization/backends/lapack/thomas.cpp
/// @brief LAPACK tridiagonal solver via LAPACKE_dgtsv.
#include "../seq/impl.hpp"
#include "impl.hpp"
#include <stdexcept>
#include <string>
#include <vector>
 
#if defined(NUMERICS_HAS_LAPACK)
    #include <lapacke.h>
#endif
 
namespace num::backends::lapack {
 
void thomas(const Vector& a,
            const Vector& b,
            const Vector& c,
            const Vector& d,
            Vector& x) {
#if defined(NUMERICS_HAS_LAPACK)
    const idx n = b.size();
    // dgtsv overwrites its inputs; work on copies
    std::vector<double> dl(a.data(), a.data() + (n - 1));
    std::vector<double> diag(b.data(), b.data() + n);
    std::vector<double> du(c.data(), c.data() + (n - 1));
    x = d;
    int info = LAPACKE_dgtsv(LAPACK_ROW_MAJOR,
                             static_cast<lapack_int>(n),
                             1,
                             dl.data(),
                             diag.data(),
                             du.data(),
                             x.data(),
                             1); // ldb = nrhs = 1 (row-major: cols of RHS)
    if (info != 0)
        throw std::runtime_error("thomas (lapack): dgtsv failed, info="
                                 + std::to_string(info));
#else
    seq::thomas(a, b, c, d, x);
#endif
}
 
} // namespace num::backends::lapack