KokkosLapack::gegqr

Defined in header: KokkosLapack_gegqr.hpp

template <class ExecutionSpace, class AMatrix, class TauArray, class InfoArray>
void gegqr(const ExecutionSpace& space, const int k, const AMatrix& A, const TauArray& Tau, const InfoArray& Info);

template <class AMatrix, class TauArray, class InfoArray>
void gegqr(const int k, const AMatrix& A, const TauArray& Tau, const InfoArray& Info);

Computes the matrix \(Q\) from the QR factorization of matrix \(A\) using the first \(k\) reflectors

\[Q=H(0)H(1)...H(k-1)\]

where \(A\) is, on input, a matrix previously factored using a call to geqrf and contains \(Q\) on output and \(Tau\) stores the associated scaling factors.

1. Overwrites \(A\) with the \(Q\) using the resources of space.

2. Same as 1. but uses the resources of the default execution space from AMatrix::execution_space.

The function will throw a runtime exception if \(k < 0\) or \(n < k\) where \(n\) is the number of columns in \(A\).

Parameters

space:

execution space instance.

k:

The number of reflectors to use when computing \(Q\).

A:

On input, matrix that contains the \(QR\) factors from a previous call to geqrf, on output, it contains the first \(k\) columns of \(Q\).

Tau:

rank-1 view of size min(M,N) that contains the scaling factors of the elementary reflectors.

Info:

rank-1 view of integers and of size 1: Info[0] = 0: successful exit; Info[0] < 0: if equal to ‘-i’, the i-th argument had an illegal value.

Type Requirements

• ExecutionSpace must be a Kokkos execution space

• AMatrix must be a Kokkos View of rank 2 that satisfies

  • Kokkos::SpaceAccessibility<ExecutionSpace, typename AMatrix::memory_space>::accessible

• Tau must be a Kokkos View of rank 1 that satisfies

  • Kokkos::SpaceAccessibility<ExecutionSpace, typename Tau::memory_space>::accessible

• Info must be a Kokkos View of rank 1 that satisfies

  • Kokkos::SpaceAccessibility<ExecutionSpace, typename Info::memory_space>::accessible

Example

// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
// SPDX-FileCopyrightText: Copyright Contributors to the Kokkos project

#include <Kokkos_Core.hpp>
#include <KokkosLapack_geqrf.hpp>
#include <KokkosLapack_gegqr.hpp>

int main(void) {
  bool correct = true;
  Kokkos::initialize();
  {
    using execution_space = Kokkos::DefaultExecutionSpace;
    using KAT             = KokkosKernels::ArithTraits<double>;

    Kokkos::View<double**, Kokkos::LayoutLeft> A("A", 3, 3), R("R", 3, 3);
    Kokkos::View<double*, Kokkos::LayoutLeft> tau("tau", 3);
    Kokkos::View<int*, Kokkos::LayoutLeft> Info("Info", 1);

    auto h_A = Kokkos::create_mirror_view(A);

    h_A(0, 0) = 12;
    h_A(0, 1) = -51;
    h_A(0, 2) = 4;
    h_A(1, 0) = 6;
    h_A(1, 1) = 167;
    h_A(1, 2) = -68;
    h_A(2, 0) = -4;
    h_A(2, 1) = 24;
    h_A(2, 2) = -41;
    Kokkos::deep_copy(A, h_A);
    Kokkos::deep_copy(R, A);

    execution_space space{};
    KokkosLapack::geqrf(space, A, tau, Info);
    KokkosLapack::gegqr(space, 3, A, tau, Info);
    Kokkos::fence();
    Kokkos::deep_copy(h_A, A);

    Kokkos::View<double**, Kokkos::DefaultHostExecutionSpace> analytic("result", 3, 3);
    analytic(0, 0) = -6. / 7;
    analytic(0, 1) = 69. / 175;
    analytic(0, 2) = 58. / 175;
    analytic(1, 0) = -3. / 7;
    analytic(1, 1) = -158. / 175;
    analytic(1, 2) = -6. / 175;
    analytic(2, 0) = 2. / 7;
    analytic(2, 1) = -6. / 35;
    analytic(2, 2) = 33. / 35;

    for (int rowIdx = 0; rowIdx < h_A.extent_int(0); ++rowIdx) {
      for (int colIdx = 0; colIdx < h_A.extent_int(1); ++colIdx) {
        if (KAT::abs(h_A(rowIdx, colIdx) - analytic(rowIdx, colIdx)) >
            10 * KAT::abs(analytic(rowIdx, colIdx)) * KAT::epsilon()) {
          std::cout << "h_A(" << rowIdx << ", " << colIdx << ")=" << h_A(rowIdx, colIdx) << ", analytic is "
                    << analytic(rowIdx, colIdx) << std::endl;
          correct = false;
        }
      }
    }
  }
  Kokkos::finalize();

  if (!correct) throw std::runtime_error("KokkosLapack::geqrf is incorrect!");
  return 0;
}