qr.cpp

Go to the documentation of this file.

/// @file linalg/factorization/backends/lapack/qr.cpp
/// @brief LAPACK QR factorization via LAPACKE_dgeqrf + LAPACKE_dorgqr.
#include "../seq/impl.hpp"
#include "impl.hpp"
#include <algorithm>
#include <stdexcept>
#include <string>
#include <vector>
 
#if defined(NUMERICS_HAS_LAPACK)
#include <lapacke.h>
#endif
 
namespace num::backends::lapack {
 
QRResult qr(const Matrix &A) {
#if defined(NUMERICS_HAS_LAPACK)
  const idx m = A.rows(), n = A.cols();
  const idx k = std::min(m, n);
 
  Matrix R = A;
  std::vector<double> tau(k);
 
  int info =
      LAPACKE_dgeqrf(LAPACK_ROW_MAJOR, static_cast<lapack_int>(m),
                     static_cast<lapack_int>(n), R.data(),
                     static_cast<lapack_int>(n), // lda = cols (row-major)
                     tau.data());
  if (info != 0)
    throw std::runtime_error("qr (lapack): dgeqrf failed, info=" +
                             std::to_string(info));
 
  // Extract upper triangle as R
  Matrix Rmat = R;
  for (idx i = 1; i < m; ++i)
    for (idx j = 0; j < std::min(i, n); ++j)
      Rmat(i, j) = 0.0;
 
  // Build full Q from stored reflectors via dorgqr
  Matrix Q(m, m, 0.0);
  for (idx j = 0; j < k; ++j)
    for (idx i = 0; i < m; ++i)
      Q(i, j) = R(i, j);
 
  info = LAPACKE_dorgqr(LAPACK_ROW_MAJOR, static_cast<lapack_int>(m),
                        static_cast<lapack_int>(m), static_cast<lapack_int>(k),
                        Q.data(),
                        static_cast<lapack_int>(m), // lda = cols = m (square)
                        tau.data());
  if (info != 0)
    throw std::runtime_error("qr (lapack): dorgqr failed, info=" +
                             std::to_string(info));
 
  return {std::move(Q), std::move(Rmat)};
#else
  return seq::qr(A);
#endif
}
 
} // namespace num::backends::lapack