lu.cpp

Go to the documentation of this file.

/// @file linalg/factorization/backends/lapack/lu.cpp
/// @brief LAPACK LU factorization via LAPACKE_dgetrf.
///
/// Falls back to seq::lu() with a warning when NUMERICS_HAS_LAPACK is not set;
/// this keeps the build working but reminds the caller that LAPACK is absent.
#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 {
 
LUResult lu(const Matrix& A) {
#if defined(NUMERICS_HAS_LAPACK)
    const idx n = A.rows();
    LUResult f;
    f.LU = A;
    f.piv.resize(n);
    f.singular = false;
 
    std::vector<lapack_int> ipiv(n);
    int info = LAPACKE_dgetrf(LAPACK_ROW_MAJOR,
                              static_cast<lapack_int>(n),
                              static_cast<lapack_int>(n),
                              f.LU.data(),
                              static_cast<lapack_int>(n),
                              ipiv.data());
    if (info < 0)
        throw std::runtime_error("lu (lapack): dgetrf argument error, info="
                                 + std::to_string(info));
    if (info > 0)
        f.singular = true;
 
    // Convert 1-based LAPACK IPIV to 0-based sequential swap targets
    for (idx k = 0; k < n; ++k)
        f.piv[k] = static_cast<idx>(ipiv[k] - 1);
 
    return f;
#else
    return seq::lu(A);
#endif
}
 
} // namespace num::backends::lapack