KokkosBatched::Spr

Defined in header: KokkosBatched_Spr.hpp

template <typename ArgUplo, typename ArgTrans>
struct SerialSpr {
  template <typename ScalarType, typename XViewType, typename APViewType>
  KOKKOS_INLINE_FUNCTION static int invoke(const ScalarType alpha, const XViewType &x, const APViewType &ap);
};

template <typename MemberType, typename ArgUplo, typename ArgTrans>
struct TeamSpr {
  template <typename ScalarType, typename XViewType, typename APViewType>
  KOKKOS_INLINE_FUNCTION static int invoke(const MemberType &member, const ScalarType alpha, const XViewType &x,
                                           const APViewType &ap);
};

template <typename MemberType, typename ArgUplo, typename ArgTrans>
struct TeamVectorSpr {
  template <typename ScalarType, typename XViewType, typename APViewType>
  KOKKOS_INLINE_FUNCTION static int invoke(const MemberType &member, const ScalarType alpha, const XViewType &x,
                                           const APViewType &ap);
};

Performs one of the symmetric or hermitian rank-1 operation of matrix \(A\) by vector \(x\) with scaling factor \(alpha\), where \(A\) is stored in packed format.

\[\begin{split}\begin{align} A &= A + \alpha (x * x^T) \: \text{(if ArgTrans == KokkosBatched::Trans::Transpose)} \\ A &= A + \alpha (x * x^H) \: \text{(if ArgTrans == KokkosBatched::Trans::ConjTranspose)} \end{align}\end{split}\]

  1. If ArgTrans == KokkosBatched::Trans::Transpose, this operation is equivalent to the BLAS routine SSPR (CSPR) or DSPR (ZSPR) for single or double precision for real (complex) symmetric matrix.

  2. If ArgTrans == KokkosBatched::Trans::ConjTranspose, this operation is equivalent to the BLAS routine CHPR or ZHPR for single or double precision for complex hermitian matrix.

Parameters

alpha:

Scaling factor

x:

On input, \(x\) is a length \(n\) vector.

ap:

With ArgUplo == KokkosBatched::Uplo::Upper, the upper triangular part of \(A\) is referenced and updated. With ArgUplo == KokkosBatched::Uplo::Lower, the lower triangular part of \(A\) is referenced and updated. The matrix \(A\) is supplied in the packed format as a length \(n*(n+1)/2\) vector, where \(n\) is the size of the matrix.

Type Requirements

  • MemberType must be a Kokkos team member handle (only for TeamSpr and TeamVectorSpr)

  • ArgUplo must be one of the following:

    • KokkosBatched::Uplo::Upper to update the upper triangular part of \(A\)

    • KokkosBatched::Uplo::Lower to update the lower triangular part of \(A\)

  • ArgTrans must be one of the following:

    • KokkosBatched::Trans::Transpose to perform \(A = A + \alpha (x * x^T)\)

    • KokkosBatched::Trans::ConjTranspose to perform \(A = A + \alpha (x * x^H)\)

  • ScalarType must be a built-in floating point type (float, double, Kokkos::complex<float>, Kokkos::complex<double>)

  • XViewType must be a Kokkos View of rank 1 containing a vector \(X\)

  • APViewType must be a Kokkos View of rank 1 containing a matrix \(A\) in packed format that satisfies std::is_same_v<typename APViewType::value_type, typename APViewType::non_const_value_type>

Example

  1#include <Kokkos_Core.hpp>
  2#include <Kokkos_Random.hpp>
  3#include <KokkosBatched_Spr.hpp>
  4
  5using ExecutionSpace = Kokkos::DefaultExecutionSpace;
  6
  7/// \brief Example of batched spr
  8/// Perform a symmetric rank-1 update of matrix A by vector x, where
  9///   A: [[1, -3, -2,  0],
 10///       [0,  3, -1, -2],
 11///       [0,  0,  9,  5],
 12///       [0,  0,  0, 27]]
 13///   A_packed: [1, -3, -2, 0, 3, -1, -2, 9, 5, 27]
 14///   x: [1, 2, 3, 4]
 15///   alpha: 1.5
 16///
 17/// After, A = A + alpha * (x * x^T), it will give
 18///   A: [2.5,  0.,  9.,  2.5,  8.,  22.5,  6.,  10.,  23.,  51.]
 19///
 20int main(int /*argc*/, char** /*argv*/) {
 21  Kokkos::initialize();
 22  {
 23    using View2DType = Kokkos::View<double**, ExecutionSpace>;
 24    const int Nb = 10, n = 4;
 25    const int packed_size = n * (n + 1) / 2;
 26
 27    // Matrix A in packed storage
 28    View2DType ap("ap", Nb, packed_size), Ref("Ref", Nb, packed_size);
 29
 30    // Vector x
 31    View2DType x("x", Nb, n);
 32
 33    // Initialize A and x
 34    auto h_ap  = Kokkos::create_mirror_view(ap);
 35    auto h_x   = Kokkos::create_mirror_view(x);
 36    auto h_Ref = Kokkos::create_mirror_view(Ref);
 37
 38    // Upper triangular matrix
 39    for (int ib = 0; ib < Nb; ib++) {
 40      // Fill vector x
 41      for (int i = 0; i < n; i++) {
 42        h_x(ib, i) = i + 1;
 43      }
 44
 45      // Fill the matrix A
 46      h_ap(ib, 0) = 1;
 47      h_ap(ib, 1) = -3;
 48      h_ap(ib, 2) = 3;
 49      h_ap(ib, 3) = -2;
 50      h_ap(ib, 4) = -1;
 51      h_ap(ib, 5) = 9;
 52      h_ap(ib, 6) = 0;
 53      h_ap(ib, 7) = -2;
 54      h_ap(ib, 8) = 5;
 55      h_ap(ib, 9) = 27;
 56
 57      h_Ref(ib, 0) = 2.5;
 58      h_Ref(ib, 1) = 0;
 59      h_Ref(ib, 2) = 9;
 60      h_Ref(ib, 3) = 2.5;
 61      h_Ref(ib, 4) = 8;
 62      h_Ref(ib, 5) = 22.5;
 63      h_Ref(ib, 6) = 6;
 64      h_Ref(ib, 7) = 10;
 65      h_Ref(ib, 8) = 23;
 66      h_Ref(ib, 9) = 51;
 67    }
 68    Kokkos::deep_copy(ap, h_ap);
 69    Kokkos::deep_copy(x, h_x);
 70
 71    // Compute A = A + alpha * (x * x^T)
 72    const double alpha = 1.5;
 73    ExecutionSpace exec;
 74    using policy_type = Kokkos::RangePolicy<ExecutionSpace, Kokkos::IndexType<int>>;
 75    policy_type policy{exec, 0, Nb};
 76    Kokkos::parallel_for(
 77        "syr", policy, KOKKOS_LAMBDA(int ib) {
 78          auto sub_ap = Kokkos::subview(ap, ib, Kokkos::ALL);
 79          auto sub_x  = Kokkos::subview(x, ib, Kokkos::ALL);
 80
 81          // A = A + alpha * (x * x^T)
 82          KokkosBatched::SerialSpr<KokkosBatched::Uplo::Upper, KokkosBatched::Trans::Transpose>::invoke(alpha, sub_x,
 83                                                                                                        sub_ap);
 84        });
 85
 86    // Confirm that the results are correct
 87    Kokkos::deep_copy(h_ap, ap);
 88    bool correct = true;
 89    double eps   = 1.0e-12;
 90    for (int ib = 0; ib < Nb; ib++) {
 91      for (int i = 0; i < packed_size; i++) {
 92        if (Kokkos::abs(h_ap(ib, i) - h_Ref(ib, i)) > eps) correct = false;
 93      }
 94    }
 95
 96    if (correct) {
 97      std::cout << "spr works correctly!" << std::endl;
 98    }
 99  }
100  Kokkos::finalize();
101}

output:

spr works correctly!