28 Numerics library [numerics]

28.9 Basic linear algebra algorithms [linalg]

28.9.14 BLAS 2 algorithms [linalg.algs.blas2]

28.9.14.1 General matrix-vector product [linalg.algs.blas2.gemv]

[Note 1:

These functions correspond to the BLAS function xGEMV.

— end note]

The following elements apply to all functions in [linalg.algs.blas2.gemv].

Mandates:

(3.1)
possibly-multipliable<decltype(A), decltype(x), decltype(y)>() is true, and
(3.2)
possibly-addable<decltype(x), decltype(y), decltype(z)>() is true for those overloads that take a z parameter.

Preconditions:

(4.1)
multipliable(A,x,y) is true, and
(4.2)
addable(x,y,z) is true for those overloads that take a z parameter.

Complexity:

O (x.extent(0) \times A.extent(1))

🔗

template<in-matrix InMat, in-vector InVec, out-vector OutVec>
  void matrix_vector_product(InMat A, InVec x, OutVec y);
template<class ExecutionPolicy, in-matrix InMat, in-vector InVec, out-vector OutVec>
  void matrix_vector_product(ExecutionPolicy&& exec, InMat A, InVec x, OutVec y);

These functions perform an overwriting matrix-vector product.

Effects: Computes

y = A x

[Example 1: constexpr size_t num_rows = 5; constexpr size_t num_cols = 6; // y = 3.0 * A * x void scaled_matvec_1(mdspan<double, extents<size_t, num_rows, num_cols>> A, mdspan<double, extents<size_t, num_cols>> x, mdspan<double, extents<size_t, num_rows>> y) { matrix_vector_product(scaled(3.0, A), x, y); } // z = 7.0 times the transpose of A, times y void scaled_transposed_matvec(mdspan<double, extents<size_t, num_rows, num_cols>> A, mdspan<double, extents<size_t, num_rows>> y, mdspan<double, extents<size_t, num_cols>> z) { matrix_vector_product(scaled(7.0, transposed(A)), y, z); } — end example]

🔗

  template<in-matrix InMat, in-vector InVec1, in-vector InVec2, out-vector OutVec>
    void matrix_vector_product(InMat A, InVec1 x, InVec2 y, OutVec z);
  template<class ExecutionPolicy,
           in-matrix InMat, in-vector InVec1, in-vector InVec2, out-vector OutVec>
    void matrix_vector_product(ExecutionPolicy&& exec,
                               InMat A, InVec1 x, InVec2 y, OutVec z);

These functions performs an updating matrix-vector product.

Effects: Computes

z = y + A x

Remarks: z may alias y.

[Example 2: // y = 3.0 * A * x + 2.0 * y void scaled_matvec_2(mdspan<double, extents<size_t, num_rows, num_cols>> A, mdspan<double, extents<size_t, num_cols>> x, mdspan<double, extents<size_t, num_rows>> y) { matrix_vector_product(scaled(3.0, A), x, scaled(2.0, y), y); } — end example]

28.9.14.2 Symmetric matrix-vector product [linalg.algs.blas2.symv]

[Note 1:

These functions correspond to the BLAS functions xSYMV and xSPMV[bib].

— end note]

The following elements apply to all functions in [linalg.algs.blas2.symv].

Mandates:

(3.1)
If InMat has layout_blas_packed layout, then the layout's Triangle template argument has the same type as the function's Triangle template argument;
(3.2)
compatible-static-extents<decltype(A), decltype(A)>(0, 1) is true;
(3.3)
possibly-multipliable<decltype(A), decltype(x), decltype(y)>() is true; and
(3.4)
possibly-addable<decltype(x), decltype(y), decltype(z)>() is true for those overloads that take a z parameter.

Preconditions:

(4.1)
A.extent(0) equals A.extent(1),
(4.2)
multipliable(A,x,y) is true, and
(4.3)
addable(x,y,z) is true for those overloads that take a z parameter.

Complexity:

O (x.extent(0) \times A.extent(1))

🔗

template<in-matrix InMat, class Triangle, in-vector InVec, out-vector OutVec>
  void symmetric_matrix_vector_product(InMat A, Triangle t, InVec x, OutVec y);
template<class ExecutionPolicy,
         in-matrix InMat, class Triangle, in-vector InVec, out-vector OutVec>
  void symmetric_matrix_vector_product(ExecutionPolicy&& exec,
                                       InMat A, Triangle t, InVec x, OutVec y);

These functions perform an overwriting symmetric matrix-vector product, taking into account the Triangle parameter that applies to the symmetric matrix A ([linalg.general]).

Effects: Computes

y = A x

🔗

template<in-matrix InMat, class Triangle, in-vector InVec1, in-vector InVec2, out-vector OutVec>
  void symmetric_matrix_vector_product(InMat A, Triangle t, InVec1 x, InVec2 y, OutVec z);
template<class ExecutionPolicy,
         in-matrix InMat, class Triangle, in-vector InVec1, in-vector InVec2, out-vector OutVec>
  void symmetric_matrix_vector_product(ExecutionPolicy&& exec,
                                       InMat A, Triangle t, InVec1 x, InVec2 y, OutVec z);

These functions perform an updating symmetric matrix-vector product, taking into account the Triangle parameter that applies to the symmetric matrix A ([linalg.general]).

Effects: Computes

z = y + A x

Remarks: z may alias y.

28.9.14.3 Hermitian matrix-vector product [linalg.algs.blas2.hemv]

[Note 1:

These functions correspond to the BLAS functions xHEMV and xHPMV[bib].

— end note]

The following elements apply to all functions in [linalg.algs.blas2.hemv].

Mandates:

(3.1)
If InMat has layout_blas_packed layout, then the layout's Triangle template argument has the same type as the function's Triangle template argument;
(3.2)
compatible-static-extents<decltype(A), decltype(A)>(0, 1) is true;
(3.3)
possibly-multipliable<decltype(A), decltype(x), decltype(y)>() is true; and
(3.4)
possibly-addable<decltype(x), decltype(y), decltype(z)>() is true for those overloads that take a z parameter.

Preconditions:

(4.1)
A.extent(0) equals A.extent(1),
(4.2)
multipliable(A, x, y) is true, and
(4.3)
addable(x, y, z) is true for those overloads that take a z parameter.

Complexity:

O (x.extent(0) \times A.extent(1))

🔗

template<in-matrix InMat, class Triangle, in-vector InVec, out-vector OutVec>
  void hermitian_matrix_vector_product(InMat A, Triangle t, InVec x, OutVec y);
template<class ExecutionPolicy,
         in-matrix InMat, class Triangle, in-vector InVec, out-vector OutVec>
  void hermitian_matrix_vector_product(ExecutionPolicy&& exec,
                                       InMat A, Triangle t, InVec x, OutVec y);

These functions perform an overwriting Hermitian matrix-vector product, taking into account the Triangle parameter that applies to the Hermitian matrix A ([linalg.general]).

Effects: Computes

y = A x

🔗

template<in-matrix InMat, class Triangle, in-vector InVec1, in-vector InVec2, out-vector OutVec>
  void hermitian_matrix_vector_product(InMat A, Triangle t, InVec1 x, InVec2 y, OutVec z);
template<class ExecutionPolicy,
         in-matrix InMat, class Triangle, in-vector InVec1, in-vector InVec2, out-vector OutVec>
  void hermitian_matrix_vector_product(ExecutionPolicy&& exec,
                                       InMat A, Triangle t, InVec1 x, InVec2 y, OutVec z);

These functions perform an updating Hermitian matrix-vector product, taking into account the Triangle parameter that applies to the Hermitian matrix A ([linalg.general]).

Effects: Computes

z = y + A x

Remarks: z may alias y.

28.9.14.4 Triangular matrix-vector product [linalg.algs.blas2.trmv]

[Note 1:

These functions correspond to the BLAS functions xTRMV and xTPMV[bib].

— end note]

The following elements apply to all functions in [linalg.algs.blas2.trmv].

Mandates:

(3.1)
If InMat has layout_blas_packed layout, then the layout's Triangle template argument has the same type as the function's Triangle template argument;
(3.2)
compatible-static-extents<decltype(A), decltype(A)>(0, 1) is true;
(3.3)
compatible-static-extents<decltype(A), decltype(y)>(0, 0) is true;
(3.4)
compatible-static-extents<decltype(A), decltype(x)>(0, 0) is true for those overloads that take an x parameter; and
(3.5)
compatible-static-extents<decltype(A), decltype(z)>(0, 0) is true for those overloads that take a z parameter.

Preconditions:

(4.1)
A.extent(0) equals A.extent(1),
(4.2)
A.extent(0) equals y.extent(0),
(4.3)
A.extent(0) equals x.extent(0) for those overloads that take an x parameter, and
(4.4)
A.extent(0) equals z.extent(0) for those overloads that take a z parameter.

🔗

template<in-matrix InMat, class Triangle, class DiagonalStorage, in-vector InVec,
         out-vector OutVec>
  void triangular_matrix_vector_product(InMat A, Triangle t, DiagonalStorage d, InVec x, OutVec y);
template<class ExecutionPolicy,
         in-matrix InMat, class Triangle, class DiagonalStorage, in-vector InVec,
         out-vector OutVec>
  void triangular_matrix_vector_product(ExecutionPolicy&& exec,
                                        InMat A, Triangle t, DiagonalStorage d, InVec x, OutVec y);

These functions perform an overwriting triangular matrix-vector product, taking into account the Triangle and DiagonalStorage parameters that apply to the triangular matrix A ([linalg.general]).

Effects: Computes

y = A x

Complexity:

O (x.extent(0) \times A.extent(1))

🔗

template<in-matrix InMat, class Triangle, class DiagonalStorage, inout-vector InOutVec>
  void triangular_matrix_vector_product(InMat A, Triangle t, DiagonalStorage d, InOutVec y);
template<class ExecutionPolicy,
         in-matrix InMat, class Triangle, class DiagonalStorage, inout-vector InOutVec>
  void triangular_matrix_vector_product(ExecutionPolicy&& exec,
                                        InMat A, Triangle t, DiagonalStorage d, InOutVec y);

These functions perform an in-place triangular matrix-vector product, taking into account the Triangle and DiagonalStorage parameters that apply to the triangular matrix A ([linalg.general]).

[Note 2:

Performing this operation in place hinders parallelization.

However, other ExecutionPolicy specific optimizations, such as vectorization, are still possible.

— end note]

Effects: Computes a vector

y^{'}

such that

y^{'} = A y

, and assigns each element of

y^{'}

to the corresponding element of y.

Complexity:

O (y.extent(0) \times A.extent(1))

🔗

template<in-matrix InMat, class Triangle, class DiagonalStorage,
         in-vector InVec1, in-vector InVec2, out-vector OutVec>
  void triangular_matrix_vector_product(InMat A, Triangle t, DiagonalStorage d,
                                        InVec1 x, InVec2 y, OutVec z);
template<class ExecutionPolicy, in-matrix InMat, class Triangle, class DiagonalStorage,
         in-vector InVec1, in-vector InVec2, out-vector OutVec>
  void triangular_matrix_vector_product(ExecutionPolicy&& exec,
                                        InMat A, Triangle t, DiagonalStorage d,
                                        InVec1 x, InVec2 y, OutVec z);

These functions perform an updating triangular matrix-vector product, taking into account the Triangle and DiagonalStorage parameters that apply to the triangular matrix A ([linalg.general]).

Effects: Computes

z = y + A x

Complexity:

O (x.extent(0) \times A.extent(1))

Remarks: z may alias y.

28.9.14.5 Solve a triangular linear system [linalg.algs.blas2.trsv]

[Note 1:

These functions correspond to the BLAS functions xTRSV and xTPSV[bib].

— end note]

The following elements apply to all functions in [linalg.algs.blas2.trsv].

Mandates:

(3.1)
If InMat has layout_blas_packed layout, then the layout's Triangle template argument has the same type as the function's Triangle template argument;
(3.2)
compatible-static-extents<decltype(A), decltype(A)>(0, 1) is true;
(3.3)
compatible-static-extents<decltype(A), decltype(b)>(0, 0) is true; and
(3.4)
compatible-static-extents<decltype(A), decltype(x)>(0, 0) is true for those overloads that take an x parameter.

Preconditions:

(4.1)
A.extent(0) equals A.extent(1),
(4.2)
A.extent(0) equals b.extent(0), and
(4.3)
A.extent(0) equals x.extent(0) for those overloads that take an x parameter.

🔗

  template<in-matrix InMat, class Triangle, class DiagonalStorage,
           in-vector InVec, out-vector OutVec, class BinaryDivideOp>
    void triangular_matrix_vector_solve(InMat A, Triangle t, DiagonalStorage d,
                                        InVec b, OutVec x, BinaryDivideOp divide);
  template<class ExecutionPolicy, in-matrix InMat, class Triangle, class DiagonalStorage,
           in-vector InVec, out-vector OutVec, class BinaryDivideOp>
    void triangular_matrix_vector_solve(ExecutionPolicy&& exec,
                                        InMat A, Triangle t, DiagonalStorage d,
                                        InVec b, OutVec x, BinaryDivideOp divide);

These functions perform a triangular solve, taking into account the Triangle and DiagonalStorage parameters that apply to the triangular matrix A ([linalg.general]).

Effects: Computes a vector

x^{'}

such that

b = A x^{'}

, and assigns each element of

x^{'}

to the corresponding element of x.

If no such

x^{'}

exists, then the elements of x are valid but unspecified.

Complexity:

O (A.extent(1) \times b.extent(0))

🔗

template<in-matrix InMat, class Triangle, class DiagonalStorage,
         in-vector InVec, out-vector OutVec>
  void triangular_matrix_vector_solve(InMat A, Triangle t, DiagonalStorage d, InVec b, OutVec x);

Effects: Equivalent to: triangular_matrix_vector_solve(A, t, d, b, x, divides<void>{});

🔗

template<class ExecutionPolicy, in-matrix InMat, class Triangle, class DiagonalStorage,
         in-vector InVec, out-vector OutVec>
  void triangular_matrix_vector_solve(ExecutionPolicy&& exec,
                                      InMat A, Triangle t, DiagonalStorage d, InVec b, OutVec x);

Effects: Equivalent to: triangular_matrix_vector_solve(std::forward<ExecutionPolicy>(exec), A, t, d, b, x, divides<void>{});

🔗

template<in-matrix InMat, class Triangle, class DiagonalStorage,
         inout-vector InOutVec, class BinaryDivideOp>
  void triangular_matrix_vector_solve(InMat A, Triangle t, DiagonalStorage d,
                                      InOutVec b, BinaryDivideOp divide);
template<class ExecutionPolicy, in-matrix InMat, class Triangle, class DiagonalStorage,
         inout-vector InOutVec, class BinaryDivideOp>
  void triangular_matrix_vector_solve(ExecutionPolicy&& exec,
                                      InMat A, Triangle t, DiagonalStorage d,
                                      InOutVec b, BinaryDivideOp divide);

These functions perform an in-place triangular solve, taking into account the Triangle and DiagonalStorage parameters that apply to the triangular matrix A ([linalg.general]).

[Note 2:

Performing triangular solve in place hinders parallelization.

However, other ExecutionPolicy specific optimizations, such as vectorization, are still possible.

— end note]

Effects: Computes a vector

x^{'}

such that

b = A x^{'}

, and assigns each element of

x^{'}

to the corresponding element of b.

If no such

x^{'}

exists, then the elements of b are valid but unspecified.

Complexity:

O (A.extent(1) \times b.extent(0))

🔗

template<in-matrix InMat, class Triangle, class DiagonalStorage, inout-vector InOutVec>
  void triangular_matrix_vector_solve(InMat A, Triangle t, DiagonalStorage d, InOutVec b);

Effects: Equivalent to: triangular_matrix_vector_solve(A, t, d, b, divides<void>{});

🔗

template<class ExecutionPolicy,
         in-matrix InMat, class Triangle, class DiagonalStorage, inout-vector InOutVec>
  void triangular_matrix_vector_solve(ExecutionPolicy&& exec,
                                      InMat A, Triangle t, DiagonalStorage d, InOutVec b);

Effects: Equivalent to: triangular_matrix_vector_solve(std::forward<ExecutionPolicy>(exec), A, t, d, b, divides<void>{});

28.9.14.6 Rank-1 (outer product) update of a matrix [linalg.algs.blas2.rank1]

🔗

template<in-vector InVec1, in-vector InVec2, inout-matrix InOutMat>
  void matrix_rank_1_update(InVec1 x, InVec2 y, InOutMat A);
template<class ExecutionPolicy, in-vector InVec1, in-vector InVec2, inout-matrix InOutMat>
  void matrix_rank_1_update(ExecutionPolicy&& exec, InVec1 x, InVec2 y, InOutMat A);

These functions perform a nonsymmetric nonconjugated rank-1 update.

[Note 1:

These functions correspond to the BLAS functions xGER (for real element types) and xGERU (for complex element types)[bib].

— end note]

Mandates: possibly-multipliable<InOutMat, InVec2, InVec1>() is true.

Preconditions: multipliable(A, y, x) is true.

Effects: Computes a matrix

A^{'}

such that

A^{'} = A + x y^{T}

, and assigns each element of

A^{'}

to the corresponding element of A.

Complexity:

O (x.extent(0) \times y.extent(0))

🔗

template<in-vector InVec1, in-vector InVec2, inout-matrix InOutMat>
  void matrix_rank_1_update_c(InVec1 x, InVec2 y, InOutMat A);
template<class ExecutionPolicy, in-vector InVec1, in-vector InVec2, inout-matrix InOutMat>
  void matrix_rank_1_update_c(ExecutionPolicy&& exec, InVec1 x, InVec2 y, InOutMat A);

These functions perform a nonsymmetric conjugated rank-1 update.

[Note 2:

These functions correspond to the BLAS functions xGER (for real element types) and xGERC (for complex element types)[bib].

— end note]

Effects:

(7.1)
For the overloads without an ExecutionPolicy argument, equivalent to: matrix_rank_1_update(x, conjugated(y), A);
(7.2)
otherwise, equivalent to: matrix_rank_1_update(std::forward<ExecutionPolicy>(exec), x, conjugated(y), A);

28.9.14.7 Symmetric or Hermitian Rank-1 (outer product) update of a matrix [linalg.algs.blas2.symherrank1]

[Note 1:

These functions correspond to the BLAS functions xSYR, xSPR, xHER, and xHPR[bib].

They have overloads taking a scaling factor alpha, because it would be impossible to express the update

A = A - x x^{T}

otherwise.

— end note]

The following elements apply to all functions in [linalg.algs.blas2.symherrank1].

Mandates:

(3.1)
If InOutMat has layout_blas_packed layout, then the layout's Triangle template argument has the same type as the function's Triangle template argument;
(3.2)
compatible-static-extents<decltype(A), decltype(A)>(0, 1) is true; and
(3.3)
compatible-static-extents<decltype(A), decltype(x)>(0, 0) is true.

Preconditions:

(4.1)
A.extent(0) equals A.extent(1), and
(4.2)
A.extent(0) equals x.extent(0).

Complexity:

O (x.extent(0) \times x.extent(0))

🔗

template<class Scalar, in-vector InVec, possibly-packed-inout-matrix InOutMat, class Triangle>
  void symmetric_matrix_rank_1_update(Scalar alpha, InVec x, InOutMat A, Triangle t);
template<class ExecutionPolicy,
         class Scalar, in-vector InVec, possibly-packed-inout-matrix InOutMat, class Triangle>
  void symmetric_matrix_rank_1_update(ExecutionPolicy&& exec,
                                      Scalar alpha, InVec x, InOutMat A, Triangle t);

These functions perform a symmetric rank-1 update of the symmetric matrix A, taking into account the Triangle parameter that applies to A ([linalg.general]).

Effects: Computes a matrix

A^{'}

such that

A^{'} = A + α x x^{T}

, where the scalar α is alpha, and assigns each element of

A^{'}

to the corresponding element of A.

🔗

template<in-vector InVec, possibly-packed-inout-matrix InOutMat, class Triangle>
  void symmetric_matrix_rank_1_update(InVec x, InOutMat A, Triangle t);
template<class ExecutionPolicy,
         in-vector InVec, possibly-packed-inout-matrix InOutMat, class Triangle>
  void symmetric_matrix_rank_1_update(ExecutionPolicy&& exec, InVec x, InOutMat A, Triangle t);

These functions perform a symmetric rank-1 update of the symmetric matrix A, taking into account the Triangle parameter that applies to A ([linalg.general]).

Effects: Computes a matrix

A^{'}

such that

A^{'} = A + x x^{T}

and assigns each element of

A^{'}

to the corresponding element of A.

🔗

template<class Scalar, in-vector InVec, possibly-packed-inout-matrix InOutMat, class Triangle>
  void hermitian_matrix_rank_1_update(Scalar alpha, InVec x, InOutMat A, Triangle t);
template<class ExecutionPolicy,
         class Scalar, in-vector InVec, possibly-packed-inout-matrix InOutMat, class Triangle>
  void hermitian_matrix_rank_1_update(ExecutionPolicy&& exec,
                                      Scalar alpha, InVec x, InOutMat A, Triangle t);

These functions perform a Hermitian rank-1 update of the Hermitian matrix A, taking into account the Triangle parameter that applies to A ([linalg.general]).

Effects: Computes

A^{'}

such that

A^{'} = A + α x x^{H}

, where the scalar α is alpha, and assigns each element of

A^{'}

to the corresponding element of A.

🔗

template<in-vector InVec, possibly-packed-inout-matrix InOutMat, class Triangle>
  void hermitian_matrix_rank_1_update(InVec x, InOutMat A, Triangle t);
template<class ExecutionPolicy,
         in-vector InVec, possibly-packed-inout-matrix InOutMat, class Triangle>
  void hermitian_matrix_rank_1_update(ExecutionPolicy&& exec, InVec x, InOutMat A, Triangle t);

These functions perform a Hermitian rank-1 update of the Hermitian matrix A, taking into account the Triangle parameter that applies to A ([linalg.general]).

Effects: Computes a matrix

A^{'}

such that

A^{'} = A + x x^{H}

and assigns each element of

A^{'}

to the corresponding element of A.

28.9.14.8 Symmetric and Hermitian rank-2 matrix updates [linalg.algs.blas2.rank2]

[Note 1:

These functions correspond to the BLAS functions xSYR2,xSPR2, xHER2 and xHPR2[bib].

— end note]

The following elements apply to all functions in [linalg.algs.blas2.rank2].

Mandates:

(3.1)
If InOutMat has layout_blas_packed layout, then the layout's Triangle template argument has the same type as the function's Triangle template argument;
(3.2)
compatible-static-extents<decltype(A), decltype(A)>(0, 1) is true; and
(3.3)
possibly-multipliable<decltype(A), decltype(x), decltype(y)>() is true.

Preconditions:

(4.1)
A.extent(0) equals A.extent(1), and
(4.2)
multipliable(A, x, y) is true.

Complexity:

O (x.extent(0) \times y.extent(0))

🔗

template<in-vector InVec1, in-vector InVec2,
         possibly-packed-inout-matrix InOutMat, class Triangle>
  void symmetric_matrix_rank_2_update(InVec1 x, InVec2 y, InOutMat A, Triangle t);
template<class ExecutionPolicy, in-vector InVec1, in-vector InVec2,
         possibly-packed-inout-matrix InOutMat, class Triangle>
  void symmetric_matrix_rank_2_update(ExecutionPolicy&& exec,
                                      InVec1 x, InVec2 y, InOutMat A, Triangle t);

These functions perform a symmetric rank-2 update of the symmetric matrix A, taking into account the Triangle parameter that applies to A ([linalg.general]).

Effects: Computes

A^{'}

such that

A^{'} = A + x y^{T} + y x^{T}

and assigns each element of

A^{'}

to the corresponding element of A.

🔗

template<in-vector InVec1, in-vector InVec2,
         possibly-packed-inout-matrix InOutMat, class Triangle>
  void hermitian_matrix_rank_2_update(InVec1 x, InVec2 y, InOutMat A, Triangle t);
template<class ExecutionPolicy, in-vector InVec1, in-vector InVec2,
         possibly-packed-inout-matrix InOutMat, class Triangle>
  void hermitian_matrix_rank_2_update(ExecutionPolicy&& exec,
                                      InVec1 x, InVec2 y, InOutMat A, Triangle t);

These functions perform a Hermitian rank-2 update of the Hermitian matrix A, taking into account the Triangle parameter that applies to A ([linalg.general]).

Effects: Computes

A^{'}

such that

A^{'} = A + x y^{H} + y x^{H}

and assigns each element of

A^{'}

to the corresponding element of A.