KokkosLapack::potrs

Defined in header: KokkosLapack_potrs.hpp

// Overload 1: with explicit execution space
template <class ExecutionSpace, class AViewType, class BViewType>
void potrs(const ExecutionSpace& space, const char uplo[], const AViewType& A, BViewType& B);

// Overload 2: uses AViewType::execution_space
template <class AViewType, class BViewType>
void potrs(const char uplo[], const AViewType& A, BViewType& B);

Solves a system of linear equations \(A X = B\) where \(A\) is a symmetric (or Hermitian) positive definite matrix whose Cholesky factorization has already been computed by KokkosLapack::potrf.

Given the factorization produced by potrf:

\[A = U^H U \quad \text{if uplo = 'U'}, \qquad A = L L^H \quad \text{if uplo = 'L'},\]

potrs solves for \(X\) by two triangular solves, overwriting \(B\) with the solution \(X\).

Parameters

space:

execution space instance (overload 1 only).

uplo:

'U' if the upper triangular factor is stored in A; 'L' if the lower triangular factor is stored.

A:

On entry, the triangular Cholesky factor of \(A\) as returned by potrf — upper triangle if uplo = 'U', lower triangle if uplo = 'L'. Not modified. Dimensions: (lda, n).

B:

On entry, the right-hand side matrix of dimensions (ldb, nrhs). On exit, overwritten with the solution matrix \(X\).

Type Requirements

• ExecutionSpace must be a Kokkos execution space

• AViewType must be a Kokkos View satisfying:

  • AViewType::rank == 2,

  • AViewType::value_type is const Scalar for a supported scalar type (float, double, Kokkos::complex<float>, Kokkos::complex<double>),

  • AViewType::array_layout is Kokkos::LayoutLeft,

  • the memory space of AViewType is accessible from ExecutionSpace.

• BViewType must be a Kokkos View satisfying:

  • BViewType::rank == 2,

  • BViewType::value_type is a (non-const) scalar type matching that of AViewType,

  • BViewType::array_layout is Kokkos::LayoutLeft,

  • the memory space of BViewType is accessible from ExecutionSpace.

Example

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

#include <Kokkos_Core.hpp>
#include <KokkosLapack_potrf.hpp>
#include <KokkosLapack_potrs.hpp>
#include <KokkosKernels_TestMatrixUtils.hpp>

#include <iostream>
#include <vector>

int main(int argc, char* argv[]) {
  Kokkos::initialize(argc, argv);
  {
    // Solve A*X = B for a 4x4 symmetric positive definite (SPD) matrix A
    // and two right-hand sides.
    //
    // Strategy:
    //   1. Choose a known X_exact and compute B = A * X_exact.
    //   2. Factor A with potrf (Cholesky): A = L * L^H.
    //   3. Solve with potrs:              L * L^H * X = B  ->  X overwrites B.
    //   4. Verify X ≈ X_exact.
    //
    // A (diagonally dominant, SPD):        X_exact:    B = A * X_exact:
    //   [10  1  2  1]                        [1  0]      [22   8]
    //   [ 1 10  3  2]                        [2  1]      [38  22]
    //   [ 2  3 10  1]                        [3  2]      [42  26]
    //   [ 1  2  1 10]                        [4  3]      [48  34]

    const int N = 4, nrhs = 2;
    using ViewType = Kokkos::View<double**, Kokkos::LayoutLeft>;

    ViewType A("A", N, N);
    ViewType B("B", N, nrhs);
    // clang-format off
    std::vector<std::vector<double>> A_data = {
        {10,  1,  2,  1},
        { 1, 10,  3,  2},
        { 2,  3, 10,  1},
        { 1,  2,  1, 10}};
    std::vector<std::vector<double>> B_data = {
        {22,  8},
        {38, 22},
        {42, 26},
        {48, 34}};
    // clang-format on

    Test::fill_view_from_fixture(A, A_data);
    Test::fill_view_from_fixture(B, B_data);

    // Step 1: Cholesky factorization — A is overwritten with the lower factor L
    KokkosLapack::potrf("L", A);

    // Step 2: Triangular solve — B is overwritten with the solution X
    Kokkos::View<const double**, Kokkos::LayoutLeft> Aconst(A);
    KokkosLapack::potrs("L", Aconst, B);

    auto h_B = Kokkos::create_mirror_view(B);
    Kokkos::deep_copy(h_B, B);

    // X_exact col-0 = [1, 2, 3, 4],  col-1 = [0, 1, 2, 3]
    std::cout << "Solution X (col 0 expected [1,2,3,4], col 1 expected [0,1,2,3]):\n";
    for (int i = 0; i < N; ++i)
      std::cout << "  X(" << i << ",0) = " << h_B(i, 0) << "  X(" << i << ",1) = " << h_B(i, 1) << "\n";
  }
  Kokkos::finalize();
}